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ABSTRACT 

Fourier  transform  techniques  have  been  the  favored  methods  in  the  analysis  of 
signals  and  systems.  One  major  drawback  of  Fourier  methods  is  the  difficulty  in 
analyzing  transient  and/or  non-stationary  behavior.  Recent  advances  in  the  field  of 
wavelet  theory  show  much  promise  in  alleviating  these  problems.  This  thesis  considers 
the  realizations  of  the  wavelet  decomposition  and  reconstruction  algorithms  for  the 
discrete  case.  The  major  discussion  will  involve  both  the  one  and  two  dimensional 
transforms.  We  also  present  a  multiple-phase  development  as  a  second  and  possibly  a 
preferable  method  for  decomposing  signals. 
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I.     INTRODUCTION 

The  recent  advances  of  wavelet  theory  have  spurned  much  interest  in  many 
disciplines:  earthquake  detection,  medical  EKG's,  and  multiple  resolution  signal 
processing.  We  are  most  interested  in  the  last  discipline  listed,  multiple  resolution  signal 
processing.  Using  wavelet  theory  in  this  field  opens  several  areas  of  interest.  In  the 
one-dimensional  case,  analysis  of  short-time  duration  signals  such  as  radar  pulses,  can 
possibly  lead  to  better  discrimination  of  similar  but  distinct  radar  sources.  For  the  two- 
dimensional  case,  statistical  pattern  recognition  using  wavelets  can  aid  in  identifying 
skewed  or  rotated  images.  Lastly,  for  both  dimensions,  wavelet  theory  can  be  used  in 
data  compression,  which  has  major  applications  particularly  in  transmitting  data  on  a 
bandwidth-constrained  channel.  Traditionally,  most  means  to  date  use  Fourier  methods 
to  attack  the  problems  in  the  areas  of  interest  discussed  previously.  However,  for 
transient  and/or  non-stationary  phenomena,  Fourier  techniques  give  inadequate  results  in 
the  transform  domain,  even  when  we  modify  them  to  include  a  dimension  of  time.  Thus, 
wavelets  can  serve  as  an  alternative  to  the  conventional  Fourier  transform  techniques  for 
analyzing  such  transient  and/or  non-stationary  behavior. 

This  thesis  develops  one  and  two-dimensional  invertible  discrete  wavelet  transforms 
(DWT).  386MATLAB  and  PROMATLAB  Version  3.5  are  the  software  tools  used  for 
386  or  better  PC's,  and  for  the  UNIX-based  SUN  workstations.  We  did  not  use  Fortran 
and  C  programming  languages  to  maximize  the  portability  of  the  routines,  but  such  codes 


could  be  developed  very  easily  from  this  thesis.    These  routines  can  be  used  by  other 
researchers  for  more  in-depth  analysis  of  signals  and  systems. 

Chapter  II  will  cover  the  basic  principles  of  the  wavelet  theory,  particularly  for  the 
discrete  dyadic  case.  Although  many  volumes  are  dedicated  to  the  field,  this  cursory 
outline  of  the  theory  will  provide  the  basis  for  notation  used  throughout  the  thesis. 
Chapter  III  will  entail  some  important  concepts  such  as  causality  and  spillover  effects  that 
must  be  accounted  for  in  the  DWT.  We  expand  these  ideas  further  into  the  two- 
dimensional  case  in  Chapter  IV.  Chapter  V  proposes  another  method  for  decomposing 
the  data  at  the  expense  of  physical  memory,  the  so-called  multiple-phase  DWT 
(MP/DWT).  The  MP/DWT  possibly  approximates  the  continuous  WT,  with  the  property 
of  time  invariance.  Chapter  VI  gives  a  narrative  description  of  the  MATLAB  algorithms 
developed.    Finally,  we  present  our  conclusions  in  Chapter  VII. 


H.    BASIC  WAVELET  THEORY 

Since  wavelets  are  a  relatively  new  and  less  established  field  to  most  readers,  one 
main  question  asked  is  "what  are  they?"  We  introduce  the  topic  in  here  by  an  example. 

Consider  a  function  f(t)  that  has  discrete  levels  as  shown  in  Figure  1.  We 
decompose  f(x)  into  a  lower  resolution  level.  In  other  words,  we  desire  to  smooth  (or 
average)  out  f(x),  or  to  low  pass  filter  the  function  to  a  coarser  level.  Figure  2  depicts 
the  resulting  decomposition  function  «f(x),  where  a  is  a  scaling  factor.  For  the  dyadic 
case,  we  set  a  to  two.  If  we  compare  the  original  function  f(x)  to  the  coarser  function 
af(x),  we  note  that  the  detail  is  ad(x)=f(x)-af(x)  as  shown  in  Figure  3.  If  we  had  the 
detail  and  the  decomposed  signal,  we  can  theoretically  reconstruct  f(x)  without  any  loss 
of  information.  The  decomposition  (and  the  reconstruction)  process  can  go  to  any 
arbitrary  resolution  level  m.  Usually,  we  sample  the  signal  at  a  rate  higher  than  or  equal 
to  the  Nyquist  rate,  and  set  the  resulting  sampled  data  to  be  at  the  highest  resolution  level 
(say  level  m=0).  All  other  resolution  levels  would  be  in  the  negative  direction.  This 
process  could  be  repeated  for  the  next  lower  resolution  level  by  operating  on  the  current 
level  in  an  iterated  manner. 


A-    - 


Figure  1.    Example  Function  f(x) 


Figure  2.    Blurred  «f(x)     Figure  3.    Detail  ad(x) 


The  primary  goal  is  to  find  a  set  of  orthonormal  basis  functions  that  will  do  the 
smoothing  of  the  original  level  and  retain  the  details  for  reconstruction  in  a  more 
practical  manner.  The  low  pass  filtering  process  uses  a  function  called  the  scaling 
function  <j>(x)  while  the  detail  basis  function  uses  the  wavelet  basis  function  ^(x),  both 
of  which  are  also  orthogonal  to  each  other. 

The  described  method  is  a  coarse  level  in  the  understanding  of  the  wavelet  theory. 
We  must  now  go  to  a  higher  resolution  level  of  knowledge  to  better  comprehend  the 
following  chapters. 

A.      THE  SCALING  FUNCTION  <t>(x) 

The  region  L2(R)  depicts  a  vector  space  where  a  function  f(x)  resides  and  satisfies 
the  condition  of  having  finite  energy,  or  the  inner  product  of  f(x)  with  its  complex 
conjugate  is  finite.  Consider  a  family  of  embedded  closed  subspaces,  Vm,  such  that  the 
infinitely  largest  subspace  is  L2(R)  while  the  infinitely  smallest  subspace  {0} . 


f  f(x)f(x)dx  =  (/(*)  ,/*(*)>< 


(1) 


(lower  resolution  ...  C  V.,  C  V0  C  V+1  C  ...  higher  resolution) 


Additionally,  for  each  vector  space  Vm  in  Vm+1,  let  Wm  be  an  orthogonal  complementary 
space  to  Vm  in  Vm+1  such  that  Vm+i=Vm©Wm.  We  show  the  spaces  graphically  in 
Figure  4. 


Figure  4 .    Vector  Subspaces 

Note  that  if  given  Vm  and  Wm  to  Wm+k,  we  can  obtain  the  vector  space  Vm+k+1  as  shown 
in  Figure  5. 
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Figure   5.         Another  View 


A  function  f(x)  will  lie  in  a  vector  space  Vm  if  and  only  if  it  can  be  written  as: 

m 

/»(*)  =  E  cmn  2T  4>(2mx-n)  (2) 


n  (  Z 


where  we  adjust  the  width  of  <f>(x)  for  unit  grid  spacing  at  resolution  level  m=0.    The 
inner  product  of  (f>  with  an  n-shifted  version  of  itself  must  equal  zero 
(  <</>(x),<Hx-n)>  =0). 

In  order  to  normalize  the  scaling  function  for  the  corresponding  resolution  level  in 
the  L2(R)  space,  multiply  the  scaling  function  by  a  factor  2m/2  while  also  dilating  in  the 
spatial  domain.  The  following  two  equations  simplify  the  notation  for  the  rest  of  the 
derivations. 

4>m  =  <f>m(x)  =  2?*0(2»*)  (3) 


Kn  =  4>Jt)  =  2*M2mx-n)  =  4>m(x-Tmn) 


Let  Pm  be  the  orthogonal  projection  of  a  function  in  L2(R)  onto  the  vector  space 
Vm.  Note  that  Equation  2  is  a  harmonic  series  representation  of  the  function  projected 
onto  the  Vm  space,  with  c^  denoting  the  weighting  coefficients  for  the  basis  function. 
This  is  similar  to  the  evaluation  of  the  fourier  coefficients  for  the  fourier  series  expansion 
with  its  basis  functions. 


Thus  Equation  5  can  be  viewed  as  a  convolution  of  the  function  with  the  mirrored 
scaling  function  sampled  every  2m  seconds. 

^-m*k  U.  (6) 

The  recursion  property  for  the  scaling  function  only  depends  on  two  adjacent  resolution 
levels  m  and  m+1.    Since  Vm  C  Vm+lJ  we  have: 

<J>«  =    E  (  *»  '  K+i*  >  4W  (7) 

A: 

We  define  the  filter  coefficients  in  terms  of  the  inner  product  of  Equation  7  as 

h(k)  =  /4>0(*)  *0(2x-k)dx  =  2~~2  <  4»0  ,  (J)+1  )  (8) 

We  are  now  able  to  bring  the  previous  equations  into  a  form  known  in  the  wavelet 
field  as  the  fundamental  dilation  equation.  This  is  a  two  scale  difference  equation 
relative  to  the  two  adjacent  levels. 

4>(*)  =  E  2  h<®  4>(2*-*)  <9> 


A  finite  set  of  h(n)  coefficients  gives  very  desirable  traits  for  selecting  wavelets 
with  compact  support.  Major  research  in  the  field  has  been  in  the  study  and 
determination  of  the  family  of  scaling  functions  that  are  of  compact  support  and 
orthonormal.  With  the  scaling  function  being  orthonormal,  we  define  the  following 
properties  which  are  relative  easy  to  prove  [Ref.  1:  pp.  689-691]: 


J  <J)(X)  dx 


(10) 
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/2(^)    =    1 
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«» -  i 

2 

hik)  h(k-2n)  = 

5(») 

(11) 


(12) 


(13) 


Going  into  the  spectral  domain,  the  low  pass  properties  of  the  scaling  function  are 
obvious  from  the  next  four  equations. 


$(/)=// 


(V 


where      H(f)  =  Y,  Kk)e'j2nk/ 

2/  t 


(14) 


00         (  f^ 

$(/)=  n  h 


P=i 


[2p 


(15) 


£   A(*)  =  1      -       ff(0)  =  1 


(16) 


H(f)  |2     +     |  i/(/4)  I2     -     1 


(17) 


B.      THE  WAVELET  BASIS  FUNCTION  tf(x) 

In  a  manner  analogous  to  the  development  of  the  scaling  function,  the  wavelet  basis 
function  ^(x)  determines  the  Wm  vector  space.  In  a  similar  fashion  to  Equation  2,  we 
define  the  detail  function  as  follows: 

*J&  ■  £  «*«  2?  *<2"*-n)  <18> 

The  wavelet  basis  function  is  of  the  form  of  a  constant-Q  band  pass  filter  function, 
increasing  an  octave  in  bandwidth  as  the  frequency  increases  an  octave.  Since  Wm  is 
orthogonal  to  Vm,  the  inner  product  of  each  of  the  spaces'  basis  functions  should  equal 
zero  (  <  0(x)  ,  4>(x)  >  =  0  ).  This  relationship  is  so  intertwined  that  by  knowing  a 
valid  set  of  ^'s,  a  set  of  i/^'s  can  be  derived  for  a  particular  resolution  level. 
Consequently,  knowing  the  h(n)  coefficients  will  determine  the  appropriate  coefficients 
for  the  wavelet  basis  functions  given  as  g(n). 

We  define  the  wavelet  basis  in  a  similar  fashion  as  the  dilation  property  in  Equation 
9: 

WO  =  2  £  g(k)  0(2*-*)  (19) 

k 

The  {^(x-n)}  is  an  orthonormal  set  spanning  Wm.   The  g(n)  coefficients  share  many  of 
the  properties  of  the  h(n)  sequence  with  the  following  exceptions: 


10 


U(pc)dx    =  Y,S&)  =0 


*V)  =  G 


f 

2 


$ 


2 


(20) 


(21) 


Mallat  [Ref.  1:  p.  681]  and  Daubechies  [Ref.  2:  p.  944]  chose  the  g(n)'s  in  terms 
of  the  h(n)'s  as  follows: 


g(ri)  =      (-iy-»h(l-n) 


(22) 


C.      THE  DISCRETE  WAVELET  TRANSFORM 

Note  that  Equation  6  shows  the  lower  resolution,  or  the  approximation,  signal's 
weighting  coefficients  as  a  convolution  of  the  time  reversal  of  the  h  coefficients  with  the 
signal  f(x)  sampled  every  2m  seconds.  Using  the  recursive  properties  of  the  scaling 
function,  it  is  easily  shown  that  the  approximation  coefficients  of  the  lower  resolution 
levels  also  can  be  determined  from  approximation  coefficients  of  the  higher  level.  In 
particular: 


11 


-J  ./WE  <  *..0  *~i* 

=E  <  *-0  (./w*-,<fc 

A: 

-2*£ft<B  cM+1(2n+*) 

it 

=2^M2/z-;)cm+1(/") 


(23) 


where  h(n)  =  h(-n). 

The  detail  weighting  coefficients  can  similarly  be  determined  as 

dmn  =2^Y,8(2n-j)cm(j)  (24) 

j 

Recall  that  the  d^'s  are  the  detail  coefficients  of  the  orthogonal  projection  of  f(x)  onto 
the  Wm  vector  space.  The  entire  collection  of  {dj,  or  {c.mJ  U  {d^  :  m  >  -M}  is 
called  the  discrete  wavelet  transform. 

1.     Decomposition 

Decomposing  to  any  lower  level  requires  that  the  L2(R)  function  be  initially 
convolved  with  <{>0(-x),  and  sampled  at  a  unity  interval  to  obtain  the  coefficients,  c0(n), 
as  shown  in  Equation  2,  and  graphically  depicted  next: 
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f(x) 
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t=2un 


Figure   6 


Determining  m=0  Coefficients 


Once  the  c^  coefficients  are  known,  we  use  Equations  23  and  24  in  an 
iterated  manner  to  obtain  the  set  of  coefficients  for  any  arbitrary  level.  For  the  dyadic 
case,  this  process  is  a  convolution,  a  two-times  decimation  of  the  result,  followed  by  a 
21/2  scaling  factor. 


mn 
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Figure  7.    Decomposition  Algorithm 
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2.      Reconstruction 

Recall  that  Vm+1=Vm©  Wm.  The  c^'s  reside  in  the  Vm  space  and  the  d^'s 
lie  in  the  Wm  space.  Thus,  a  projection  of  a  function  onto  the  next  higher  level  should 
equal  the  sum  of  the  projections  onto  Vm  and  Wm.  Let  Pm  and  Qm  be  the  projection 
operators  onto  the  vector  space. 

P„Jfi  =  PJfi  ®  Qjfl 

^  (25) 

n 
n 

Matching  the  terms  in  Equation  set  25  and  doing  several  transformations  of 
variables,  we  get  the  following: 


cm+i  =v/2  E  [cJl)Hn-2l)  +dm(l)g(n-2l)]  (26) 


Equation  26  shows  that  the  reconstruction  of  the  approximation  coefficients 
of  one  level  higher  can  be  realized  by  putting  zeros  between  the  lower  resolution 
coefficients,  convolving  with  the  respective  H  or  G  filter,  adding  and  then  scaling  the 
result  by  21/2.   We  show  a  block  diagram  of  this  basic  algorithm  in  Figure  8. 
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Figure  8.    Reconstruction  Algorithm 

For  the  rest  of  the  discussion,  we  will  concentrate  primarily  on  the 
approximation  and  detail  coefficients.  It  is  the  properties  of  these  coefficients  that  are 
of  importance  in  relation  to  the  DWT  process. 

The  basic  procedure  would  be  to  decompose  the  data  array,  starting  from  the 
highest  resolution  level.  Next  calculate  the  coefficients  for  the  next  lower  level,  while 
retaining  the  details,  d^,  and  using  the  resulting  approximation  coefficients  to  calculate 
the  next  lower  level.  Repeat  this  procedure  until  reaching  the  lowest  desired  level  where 
now  both  the  approximation  and  detail  coefficients  are  retained. 
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m.    ONE  DIMENSIONAL  DWT  DEVELOPMENT 

A.      DECOMPOSITION 

1.       Causality 

Recalling  that  the  decomposition  algorithm  of  Figure  7,  a  finite  data  vector 
of  the  sampled  signal,  c^,  can  be  decomposed  to  its  lower  resolution  coefficients,  c.ln  and 
d.ln.  This  process  requires  that  the  data  array  be  convolved  with  the  respective  filters. 
These  filters  are  the  time-reversed  version  of  the  H  and  G  vector  array.  At  the  output 
of  each  convolution,  there  is  a  two  times  decimation.  The  2m  scalar  multiplication 
normalizes  energy  the  in  the  L2(R)  domain. 

Consider  a  finite  array  c^  of  length  L.  Let  c^  to  be  a  non-zero  set  only  if 
the  indices,  n,  satisfy  0  <  n  <  L-l.  Also  set  the  number  of  h  coefficients  equal  to  an 
even  number  N.  By  using  Equation  23,  the  convolution  process  can  be  depicted  as 
shifting  the  h's  by  a  factor  of  two  each  time. 


C0.0 

Co., 

Co.2 

CC3 

Co.. 

c0.5 

Co.e 

Co., 

Co., 

c 

Shift  by 

h 

-2 

h 

-1 

h 

0 

h 

1 

-►   Shift  by  two 

Figure  9 .    Shift  Factor  of  2  for  N=4 
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In  an  effort  to  make  the  transform  a  causal  system,  we  desire  to  make  all  the 
non-zero  data  coefficients  have  indices  greater  or  equal  to  zero  since  negative  indices 
would  equate  to  negative  time.  This  reconstruction  gives  two  possible  sets  of  indices  for 
the  h  filter  coefficients. 

•  {h(-N+l),...,h(-2),h(-l),h(0)} 

•  {h(-N+2),...,h(-l),h(0),h(l)} 

We  allow  the  negative  indices  for  the  filter  since  we  already  know  the  filter 
coefficients  (a-priori  data).  These  two  sets  for  the  filter  coefficients  equate  to  the  two 
possible  phases  that  can  be  used  for  the  decomposition.  For  notation  purposes,  we  will 
set  the  largest  value  of  the  indices  to  equate  to  the  type  of  phase  decomposition,  phase-0 
or  phase- 1. 

By  using  Mallat  and  Daubechies  selection  of  the  g  coefficients  in  Equation 
22,  there  is  a  problem  that  occurs.  While  the  resulting  convolution  with  the  h  mirror 
filter  gives  the  values  for  indices  that  are  greater  or  equal  to  zero,  the  output  of  the  detail 
coefficients  give  values  for  negative  indices.  In  other  words,  given  we  define  the  h's 
with  phase-0  indices,  the  resulting  g  coefficient  indices  are  all  greater  or  equal  to  one, 
as  shown  in  Equation  23. 
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In  order  to  get  both  the  c^'s  and  d^'s  to  fall  into  the  right  hand  side,  we 
must  find  another  relation  for  g  in  terms  of  h,  that  satisfies  the  properties  in  Equations 
18-21.   The  following  equation  will  serve  the  purpose: 


g(ri)  =  (-iyN+3  h(-N+3-n) 


(27) 


2.       Phase  Selection 


The  decomposition  can  proceed  with  the  choice  of  phase  to  select.  The 
selection  of  the  desired  decomposition  phase  and  size  of  the  input  vector,  L,  determines 
the  total  number  of  coefficients  at  the  lower  resolution  level. 


phase-0:      |{c°_1>B}|    = 


phase-h      \{cl-lin}\ 


N+L-l 


(28) 


N+L-l 


3.       Zero-Padding  of  the  Data  Array 

The  basic  algorithm  uses  the  dot  product  between  the  h  (and  g)  coefficient 
vector  with  a  set  of  data  coefficients,  c^'s,  that  are  obtained  by  a  sliding  window  on  the 
input  coefficients,  cm+ln  vector.  The  sliding  window  shifts  to  the  right  on  a  generated 
work  vector  to  simplify  the  programming.  This  work  vector  is  a  zero-padded  version 
of  the  input  coefficient  vector.  Depending  on  the  length  of  the  input  data  array,  the 
selection  of  the  phase,  and  the  size  of  the  filters,  the  required  zero-padding  can  be 
determined  for  the  beginning  and  the  end  of  the  sequence. 


18 


If  we  select  a  phase-0  decomposition,  we  will  generate  the  work  array  by 
initially  padding  one  zero  on  the  left.  Additionally,  if  the  size  of  the  data  array  is  even, 
we  append  another  zero  to  the  end  of  the  sequence.  For  the  phase- 1  case,  we  pad  the 
work  array  with  one  zero  at  the  end  of  the  sequence  only  if  the  original  array  is  odd. 
After  determining  all  these  conditions,  we  finally  pad  the  work  array  with  N-2  zeros  both 
in  the  beginning  and  the  end  of  the  intermediate  sequence. 


Initially  pad  with  one  zero  if  the 

Phase-0   :  the  number  of  coefficients  in  the 

array  is  EVEN. 


0    0 


~i 1 1 r 


Original  coefficient  array 


N-2 


0    0     0 


Initially  pad  with  one  zero  if  the 
rhase-1   :  the  number  of  coefficients  in  the 

array  is  ODD. 


Figure  10.   Definition  of  the  Work  Array 

In  esoteric  terms,  a  phase- 1  decomposition  for  each  level  would  be  more 
efficient.  Practically,  this  savings  in  memory  is  very  insignificant,  since  the  number  of 
practical  resolution  levels  is  approximately  ["l°g2(  I  con  I  )1  •  Additionally,  this  reduction 
in  the  size  of  the  workspace  is  balanced  by  a  larger  workspace  during  the  reconstruction 
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process.  Most  importantly,  desiring  only  to  decompose  the  data  with  only  the  phase- 1 
case  for  each  level  may  miss  some  properties  of  a  signal,  which  will  be  further  discussed 
in  Chapter  V. 

4.       Energy  Determination 

One  quick  method  in  assessing  if  the  decomposition  routine  is  working 
without  doing  the  reconstruction  is  to  check  the  energy  of  the  signal  in  each  resolution 
level.  Recall  that  Equations  2  and  18  are  harmonic  series  representations.  By  using 
Poisson's  Summation  Formula  (PSF),  the  energy  in  the  function  in  the  spatial  domain  is 
equal  to  the  total  energy  stored  in  the  weighting  coefficients  of  the  transform  domain. 
This  property  holds  true  when  the  basis  functions  are  orthonormal  in  the  L2(R)  space. 
Thus,  the  energy  of  a  higher  resolution  level's  approximation  coefficients  is  equivalent 
to  the  sum  of  the  energies  of  the  approximation  and  detail  coefficients  one  resolution 
level  down. 


/  >  Cm  +  l,n    "    /  > 


cli+d2ml 


(29) 


By  normalizing  the  energy  to  the  input  resolution  level,  m=0,  a  quick  check 
of  all  the  levels  will  add  confidence  to  the  program  accuracy.  In  addition,  this  energy 
check  shows  the  numeric  truncation  errors  inherent  on  the  processing  unit  used.  The 
display  of  the  energy  will  provide  further  insight  on  possible  compression  schemes  for 
classes  of  signals,  since  the  energy  may  be  only  concentrated  at  certain  resolution  levels. 
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Now  as  an  example,  consider  a  one-dimensional  transient  signal*,  with  a 
phase- 1  decomposition  for  all  the  levels  and  using  a  Daubechies  2-tap  filter  (N=4) 
[Ref.  2:  p.  980].   The  following  diagrams  show  the  original  signal  and  the  energy  plot. 
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Figure  11.   Transient  Signal 


Professor  Michael  A.  Morgan,  Chairman  of  the  Department  of  Electrical  and 
Computer  Engineering  at  the  Naval  Postgraduate  School,  provided  the  transient  data 
included  in  the  MATLAB  routines  called  "et90.m  ."  This  signal  is  the  back-scattered 
transient  electromagnetic  field  from  a  10cm  long  thin- wire  illuminated  by  a  double- 
Gaussian  impulse,  computed  using  a  time-domain  integral  equation.  The  other  transient 
signal  provided  by  Professor  Morgan  is  "et60.m".  The  numbers  90  and  60  correspond 
to  the  angle  of  incidence  that  the  double-Gaussian  impulse  impinges  onto  the  wire's  axis. 
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Figure  12 .   Energy  in  each  Level 

The  energy  check  going  up  from  the  lowest  resolution  is  depicted  next.  Note 
that  the  Daubechies  h  filter  coefficients  are  only  significant  to  10"12  places,  thus  they 
cause  insignificant  numerical  errors  in  the  reconstruction. 


Verification  of  the  energy  in  each  resolution  level 
Level     Energy  in  Enerov  in  C  Difference 


Energy  in 
level  m-1 
c+d 


Energy  in  C 
level  m 


-7  0.000024166621151 

-6  0.000258409042490 

-5  0.001707402900115 

-4  0.015878270678385 

-3  0.197762381708774 

-2  0.474398688169256 

-1  0.936969604304463 

0  1.000000000000806 


0.000024166621151 
0.000258409042489 
0.001707402900114 
0.015878270678372 
0.197762381708634 
0.474398688168927 
0.936969604304102 
1.000000000000000 


0.000000000000000 
0.000000000000000 
0.000000000000001 
0.000000000000013 
0.000000000000140 
0.000000000000329 
0.000000000000361 
0.000000000000806 
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5.       Approximation  and  Detail  Time/Scale  Diagrams 

One  interesting  result  of  the  DWT  for  the  one  dimensional  case  is  the  display 
of  the  energy  on  the  time/scale  plane.  The  scale  can  easily  be  related  to  frequency  by 
noting  that  the  highest  resolution  is  the  sampling  frequency  at  m=0.  As  the  scale 
decreases  one  level,  we  halve  the  center  frequency  for  this  scale. 

Figure  7  shows  the  decimation  of  the  coefficients  by  a  factor  of  two.  When 
we  obtain  the  output  of  the  decomposition  process,  the  coefficient's  are  in  a  packed 
format,  meaning  that  the  appropriate  time  spacing  is  not  preserved.  Thus,  2~m-l  zeros 
must  be  padded  between  each  coefficient  for  a  particular  resolution  level  m. 

The  following  time/frequency  diagrams  show  the  transient  signal,  mentioned 
in  the  previous  section,  with  the  proper  spacing  of  the  coefficients'  energy  in  the 
transformed  domain. 

Zoon-ied      Energy     Display     of     trie     Approximation  Rr-iase:      11111111 


Figure  13. 


Mesh  Plot  for  the  Energy  of  the  Approximation 
Coefficients  to  Eight  Resolution  Levels 
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Figure  14.  Contour  Display  for  the  Energy  of  the 
Approximation  Coefficients  to  Eight  Resolution 
Levels 
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Figure  15 . 


Mesh   Plot   for   the   Energy   of   the   Detail 
Coefficients  to  Eight  Resolution  Levels 
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Figure  16. 


Contour  Display  for  the  Energy  of  the  Detail 
Coefficients  to  Eight  Resolution  Levels 


B.   RECONSTRUCTION 

Figure  8  and  Equation  25  show  the  realization  of  the  recomposition  routine  by  first 
inserting  a  zero  after  each  coefficient.  Next  we  convolve  the  array  with  their  respective 
filter,  sum,  and  finally  normalize  by  a  factor  of  2m.  This  process  is  essentially  an 
interpolation. 
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Since  we  have  defined  the  indexing  of  the  h  and  g  coefficients  during  the 
decomposition  phase,  the  recomposition  process  is  fixed  to  one  of  two  methods 
depending  on  the  decomposition  phase  used  to  get  to  that  resolution  level.  The  same 
dyadic  convolutional  algorithm  can  be  used  for  the  recomposition  after  initially  massaging 
the  input  data  arrays,  and  shifting  one  unit  to  the  right  instead  of  two.  Initially,  the  h's 
and  the  g's  are  reversed  in  order. 

Now  we  generate  the  workspace  vectors  for  the  respective  coefficients.  Append  a 
zero  AFTER  each  of  the  input  coefficient  vectors.  Then  add  two  more  zeros  at  the  end 
of  the  modified  work  vectors.  Finally,  if  the  phase- 1  decomposition  process  was  utilized 
to  obtain  the  coefficients,  append  another  zero  at  the  beginning  of  the  modified  sequence. 
The  diagram  below  summarizes  the  work  vector  definition  prior  to  running  the 
convolution  algorithm. 

N-2 


0::«  OfflOfflOfflOmOMO 


0    0 


0 


0 


Pad  one  extra  zero  for  the  phase- 1  case  ONLY 


Pad  one  zero  AFTER  every  "packed"  coefficient 


Individual  coefficients  from  the  "packed"  format 


Figure  17.   One-dimensional  Reconstruction  Work  Array 
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All  the  information  needed  for  the  reconstruction  is  the  selected  phase 
decomposition  to  get  the  c's  and  d's  for  that  resolution  level.  One  factor  to  note  is  that 
with  Equations  26,  there  is  a  possibility  that  an  extra  coefficient  be  generated  in  the 
reconstruction  process.  This  extra  coefficient  will  theoretically  equal  zero,  but 
numerically  it  may  not  and  can  be  cascaded  into  a  notable  error  if  not  accounted  for.  So 
since  we  know  the  length  of  the  input  vector  at  resolution  level  m=0,  the  size  of  the  data 
array  will  be  known  for  all  levels  since  we  retain  all  the  detail  coefficients.  Any  extra 
zero-valued  coefficient  generated  can  be  identified  as  such  and  ignored  in  going  up  to 
the  next  resolution  level. 

Errors  in  the  regeneration  will  be  totally  manifested  as  numeric  errors  of  the 
computer  processors.  The  following  diagram  shows  the  reconstructed  transient  signal 
and  absolute  error  with  the  original. 
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Figure  18.   Original  Transient  Signal  and  Reconstructed 
Versions 
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Figure  19.   Magnitude   of   the   Absolute   Error   of 
Reconstruction.   Notice  the  maximum  value. 
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IV.    TWO  DIMENSIONAL  DWT  DEVELOPMENT 

A.      INTRODUCTION 

Arguments  for  the  DWT  can  be  generated  for  any  higher  dimension.  Of  particular 
concern  is  the  two-dimensional  case  for  image  processing.  The  vector  subspaces  now 
will  reside  in  the  L2(R2)  vice  L2(R)  space.  A  two-dimensional  scaling  function  <j>(x,y) 
exists  whose  scaled  dilations  and  translations  for  a  particular  resolution  level  form  an 
orthonormal  basis  for  the  vector  subspace  Vm.  By  making  the  two-dimensional  scaling 
function  separable,  </>(x,y)=<£(x)<My),  eacn  vector  space  can  be  decomposed  as  a  tensor 
product  of  two  identical  subspaces  in  L2(R). 

V    =  V\  ®  Vl  (30) 

m  tn  m 

The  tensor  product  can  be  viewed  as  the  in-phase  and  quadrature-phase  components  in 
Vm.  The  properties  depicted  in  Figures  4  and  5  still  apply  for  the  two-dimensional  case. 

=  K©  *i)®  (wi©  K)  (31) 

=  {vL®  vl)®  (vl®  wl)®  (wl®  vl)®  (wl®  wlmi 
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Notice  that  (j)^  defines  a  set  of  orthonormal  basis  functions  in  Vm+11  and  i/^  in  Wj. 
Thus,  Equation  31  implies  that  four  sets  of  equations  form  the  orthonormal  basis  of  Vm+1 
we  summarize  in  the  following  table. 


Basis           Vector  Space         Basis                 Frequency 
Function                                Coefficient        Characteristics 

^(x)™  ^(y) 

T  m     x~*    T  m 

cm(k,l) 

Vert  Low  Pass 
Horiz  Low  Pass 

^(x)^  tanty) 

V'^W1 

T  m     x-*    *"  m 

lcL(k,l) 

Vert  Low  Pass 
Horiz  Band  Pass 

^(x)™  ^(y) 

W  l  <8>  V  ! 

2dm(k,l) 

Vert  Band  Pass 
Horiz  Low  Pass 

tWnm  ^mn(y) 

w  1  ®  w  l 

3<Uk,l) 

Vert  Band  Pass 
Horiz  Band  Pass 

Let  Pm,  lDm,  2Dm,  and  3Dm  denote  the  projection  of  a  L2(R2)  function  f(x,y),  onto 
vector  subspaces  Vm  ®  Vm,  Vm  ®  Wm,  Wm  ®  Vm,  and  Wm  <8>  Wm  respectively.  Then  in 
a  fashion  analogous  to  Equation  25,  the  vertical  and  horizontal  low  pass  operation  of 
f(x,y)  onto  resolution  level  m  is  as  follows: 


pj  *  («*o»**  Mb*)  ■  2E  E  ejM*jr***jrHn      <32> 

k,l  e  Z2 


Notice  that  the  weighting  coefficients,  cm,  are  now  indexed  in  two-dimensions. 
Similarly,  the  ldm,  2dm,  and  3dm  coefficients  can  be  determined  which  correspond  to  the 
low  horizontal  and  high  vertical,  high  horizontal  and  low  vertical,  and  high  horizontal 
and  high  vertical  frequency  details. 
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Since  the  three  two-dimensional  wavelets  and  the  one  low  pass  basis  functions  are 
separable,  the  two-dimensional  case  is  very  easily  realized.  The  weighting  coefficients 
can  be  determined  by  first  decomposing  on  the  rows,  then  the  columns  of  the  array,  or 
vice  versa.    The  general  algorithm  is  depicted  in  Equations  33-36  and  Figure  20. 


cm-i(U)  ■  2j2yEcm(ktl)hk{l-2j)hv{k-2i) 

k      I 


(33) 


k      I 


(34) 


idm.fij)  =  2Y:Y:cJk^gh(i-2j)hv(k-2i) 

k      I 


(35) 


Hn-i(U)  =  2££cw(M)sA(/-2/)sv(*-20 

k      I 


(36) 
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Figure  20.   2-D  Decomposition  Routine 
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B.      DECOMPOSITION 

We  will  use  the  convention  in  image  processing  that  the  upper  left  corner  is  the 
starting  and  the  lower  right  is  the  ending  data  point.  This  way,  we  extend  the  idea  of 
causality  into  two  dimensions,  where  the  spillover  effects  will  fall  to  the  right  and  below 
the  input  array. 

1.       Decomposition  Mask 

By  looking  at  the  'd^  coefficients'  indices  between  the  left  and  right  side  of 
the  equation,  generate  a  two-dimensional  mask  by  forming  the  outer  product  of  the 
vertical  and  horizontal  filter  coefficients. 

HhGv  -  [gf  •  [hh]  (37) 

where  [h]  =  [h(-N+2)  ...  h(0)  h(l)].  Similarly  the  masks  for  HhHv,  GhHv,  and  GhGv  also 
can  be  easily  determined. 

We  require  four  masks  to  decompose  successfully  an  image  as  suggested  in 
Figure  20  to  generate  the  coefficients  cm.,,  1dm.1, 2dm.!,  and  3dmA.  HhHv,  HhGv,  GhHv,  and 
GhGv  are  the  corresponding  filter  masks  as  shown  next: 
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Figure  21.       Mask  HhHv 
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Figure  22.       Mask  HhGv 
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Figure   23.       Mask   GhHv 
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Figure   24.       Mask  GhGv 
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These  masks  generate  the  corresponding  coefficients  by  shifting  by  a  factor 
of  two  to  the  right  and  downward  through  a  workspace  array.  With  the  one-dimensional 
case,  there  were  only  two  possible  phase  decompositions.  Now  four  possible  phases 
exist.    Obviously,  the  workspace  array  is  different  for  each  phase  decomposition. 
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Figure  25.   The  Four  Possible  2-D  Decomposition  Phases 
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2.       Zero-padding  of  the  Data  Array 

As  stated  from  the  previous  section,  the  work  space  arrays  are  different  from 
each  other  depending  on  one  of  the  four  decomposition  phases  selected,  the  size  of  the 
mask,  and  the  size  of  the  input  array.    These  considerations  must  be  noted. 

The  generation  of  the  work  array  for  the  phase-00  case  first  consists  of 
determining  if  the  input  array  rows  and/or  columns  are  odd  or  even.  If  either  of  them 
are  even,  append  a  corresponding  extra  zeros  row  or  column  to  the  bottom  or  right  of 
the  input  array,  respectively.  Given  the  number  of  rows  and  columns  for  the  mask  as 
Nv  and  Nh,  add  Nv-1  zeros  rows  above  and  Nv-2  zeros  rows  below  the  modified  input 
array.  Append  Nh-1  zeros  rows  to  the  left  and  Nh-2  zeros  rows  to  the  right  of  the 
modified  input  array. 


nv-i{ 


If  #  columns  is  even 
add  one  more  column 
of  zeros 


}Nv-2 


If  #  rows  is  even 
add  one  more  row 
of  zeros 


Figure  26.   Phase-00  Work  Array 
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For  the  phase-10  case,  if  the  original  array's  rows  are  even  and/or  if  the 
columns  are  odd,  add  one  respective  zeros  row/column  to  the  corresponding  bottom/right 
of  the  modified  array.  Next  add  Nv-1  zeros  rows  above  and  Nv-2  zeros  rows  below  the 
work  array.    Lastly,  add  Nh-2  zeros  columns  to  the  right  and  left  of  the  array. 


N„-2 


Nv-1{ 


If  #  columns  is  odd 
add  one  more  column 
of  zeros 


}Nv-2 


If  #  rows  is  even 
add  one  more  row 
of  zeros 


Figure  27.   Phase-10  Work  Array 

Alternatively  the  phase-01  case,  if  the  original  array's  rows  are  odd  and/or 
if  the  columns  are  even,  add  one  respective  zeros  row/column  to  the  corresponding 
bottom/right  of  the  modified  array.  Next  add  Nv-2  zeros  rows  above  and  below  the  work 
array.  Lastly,  add  Nh-1  zeros  columns  to  the  right  and  Nh-2  zero  columns  to  the  left  of 
the  array. 
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V2 


Nv-2{ 


If  #  columns  is  even 
add  one  more  column 
of  zeros 


}N-2 


If  #  rows  is  odd 
add  one  more  row 
of  zeros 

Figure  28.   Phase-01  Work  Array 

Lastly  for  the  phase-11  case,  if  either  of  the  input  array's  rows  or  columns 
are  odd,  then  include  a  corresponding  extra  zeros  row  or  column  to  the  bottom  or  right 
of  the  input  array,  respectively.  Given  the  number  of  rows  and  columns  for  the  mask 
as  Nv  and  Nh,  add  Nv-2  zeros  rows  above  and  below  the  modified  input  array.  Add  Nh-2 
zeros  rows  to  the  left  and  to  the  right  of  the  modified  input  array. 
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If  #  columns  is  odd 
add  one  more  column 
of  zeros 


}Nv-2 


If  #  rows  is  odd 
add  one  more  row 
of  zeros 


Figure  29.   Phase-11  Work  Array 

3.       Energy  Determination 

Similar  to  the  one-dimensional  case,  adding  up  the  squares  of  the  four  sets 
of  coefficients  at  resolution  level  m  will  equal  the  energy  in  the  coefficients  for  the  next 
higher  level,  m+ 1 .  This  is  a  natural  two-dimensional  extension  of  Poisson's  Summation 
Formula,  discussed  previously  in  Chapter  III. A. 5,  as  shown  as  follows: 


EE  lO2  ■  EE  (  K\2  ♦  l'4J2  ♦  IHJ2  ♦  IHJ2 )       (38) 

j     k  j     k 
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As  an  example,  consider  an  image  of  a  centered  square,  in  Figure  30, 
decomposed  down  to  three  resolution  levels  with  a  phase-00  HAAR  wavelet.  The  energy 
plot  in  Figure  31  shows  the  distribution  of  the  c's  and  the  four  d's  for  each  level. 
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3D  display  of  the  Image 
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Figure  30.   Square  Image 
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Figure  31.   Energy  Distribution 
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Except  for  numerical  computations  of  the  computer  processor,  all  the  energy 
from  the  lower  d  coefficients  plus  the  energy  of  the  lowest  c  coefficients  add  up  to  the 
energy  of  the  original  coefficients,  c0.  The  following  energy  check  results  for  the  two- 
dimensional  case,  normalized  to  c0.  Notice  that  there  are  no  errors  since  we  used  Haar 
h  filter  coefficients  and  that  the  two-dimensional  decomposition  uses  a  normalization 
factor  of  two  vice  the  square  root  of  two. 

Verification  of  the  energy  in  each  resolution  level 

Difference 


Level 

Energy  in 

Energy  in  C 

m 

level  m-1 
c+dl+d2+d3 

level  m 

-2      0.440942796610169       0.440942796610169     0.000000000000000 

-1      0.751059322033898       0.751059322033898     0.000000000000000 

0      1.000000000000000       1.000000000000000     0.000000000000000 


C.      RECONSTRUCTION 

Since  the  scaling  and  wavelet  basis  functions  are  separable,  the  reconstruction 
algorithm  closely  follows  the  same  process  as  in  the  one-dimensional  case.  First  each 
row/column  is  zero  padded  appropriately  and  reconstructed  with  the  one-dimensional 
technique.  Then  process  each  resulting  column/row  in  a  similar  fashion  but  with  a  one 
unit  shift  to  the  right  and  then  downward,  as  in  a  raster  scan.  We  show  this  process  in 
Figure  32. 


41 


d  (i,j 


n-l 


ml         w 

h  (i,j) 


m-l 


C  (i,j) 


•t2  — Jt2 

"*  Gh  ">  G  > 
n               v 

•t2  ->t2  - 

>H    -JG    . 

h              v 

12 

-*t2-*Gh  _ 

->  H    • 

V 

|t2 

-t2~*[V 

-►  H    > 

V 

C  (i,j) 


Figure   32.        2-D  Reconstruction 

By  taking  another  look  at  Figure  32,  a  more  efficient  method  can  be  realized.  By 
padding  each  coefficient  in  the  array  with  one  zero  in  both  the  horizontal  and  vertical 
directions,  we  can  then  slide  a  reconstruction  mask  through  the  padded  array. 

1.       Reconstruction  Mask 

We  generate  the  reconstruction  masks  in  a  similar  fashion  as  in  the 
decomposition  case,  but  with  the  important  note  that  the  coefficients  are  in  the  reverse 
order  to  accommodate  the  reconstruction  convolution  depicted  in  Figure  32  and  the 
following  equation  in  matrix  notation 


ga  -  \.KY  '  EtJ 


(39) 


where  h  =  [h(l)  h(0)  ...  h(-N+2)]. 
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Another  simpler  way  to  generate  these  masks  is  to  reverse  the  order  of  the 
decomposition  mask  in  both  the  horizontal  and  vertical  directions.  We  show  the 
comparison  of  the  two  masks  in  Figures  33  and  34  for  the  GvHh  case  only,  all  other 
reconstruction  masks  can  be  similarly  related. 
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Figure  33.   Decomposition  Mask  GhHv 
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Figure  34.   Reconstruction  Mask  HvGh 
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2.       Zero-padding  of  the  Coefficient  Arrays 

In  order  to  use  the  reconstruction  masks,  we  generate  the  work  array  for  the 
lower  level  coefficients  to  accommodate  for  the  decomposed  coefficient  array  and  mask 
size,  and  the  decomposition  phase  that  was  selected. 

The  phase-00  reconstruction  work  array  is  the  smallest  of  the  four  work 
spaces.  The  coefficient  array  is  padded  initially  with  one  zero  to  the  right  and  bottom 
of  each  of  the  individual  coefficients.  This  makes  the  modified  array  have  an  even 
number  of  rows  and  columns.  Then  append  Nh-2  columns  of  zeros  to  the  right  of  the 
array  followed  by  Nv-2  rows  of  zeros  on  the  bottom. 


o 


o 


o 


0 


o 


o 


0  0 

0  0 

0  0 

0  0 

0  0 

0  0 


N  -2  rows  of  zeros 
h 


oooooooo 

00000000 


N  -2  columns  of  zeros 
v 


Initial  Array  Coefficients  in  a  packed  format 


Q        -  Padded  zero  put  to  the  right,  bottom  or  the  bottom  right  of  the 
individual  packed  coefficients 


Figure  35.   Phase-00  Reconstruction  Work  Array 
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For  the  phase- 10  case,  add  one  extra  column  of  zeros  to  the  left  of  the  phase- 
00  work  array.  Alternatively  for  the  phase-01  case,  append  one  extra  row  of  zeros  to 
the  top  of  the  phase-00  work  space.  Lastly  for  the  phase-11  situation,  insert  both  one 
row  and  one  column  of  zeros  to  the  top  and  left  of  the  phase-00  case. 


Add  one  row  of  zeros  for 
either: 

Phase-01 

Phase-11 


Phase-00 
Reconstruction 
Work  Array 

from  the  previous 
diagram 


Add  one  column  of  zeros  for 
either: 

Phase- 10 

Phase-11 


Figure   36.        Work   Space   for   Phase-10, 10, 11   Cases 

3.       Two  Dimensional  Example 

Consider  the  diagram  of  Figure  30,  the  following  plots  are  the  resulting  mesh 
display  of  the  decomposed  coefficients  three  resolution  levels  down  (m=-3)  using  a  Haar 
wavelet  and  a  constant  phase-00  selection  for  each  level.  The  top  left  diagram  is  the 
approximation  coefficients,  cm.  Going  in  a  clockwise  manner,  the  following  plots  equate 
to  1dm,  2dm,  and  3dm  coefficients. 
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Figure   37.       Mesh  Display  of  the  Coefficients   for   level  m=-3 
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Figure   38. 


Contour   Display   of   the   Coefficients   for    level 
m=-3 
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Figures  39  and  40  compare  the  original  and  reconstructed  image  coefficients 
displayed  as  mesh  and  contour  plots. 
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Figure   39.       Reconstructed  and  Original  Mesh  Comparisons 
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Figure  40.   Reconstructed  and  Original  Contour  Comparisons 
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V.   MULTIPLE  PHASE  DWT 

For  both  the  one  and  two-dimensional  transforms,  we  selected  a  single  phase  for 
decomposing  to  the  next  lower  resolution  level.  If  we  decompose  down  to  a  fixed  level 
m=M  while  picking  different  decomposition  phases  for  each  level,  we  generate  a 
different  set  of  coefficients.  This  set  of  decomposition  coefficients  and  any  other  set  of 
decomposition  coefficients  will  successfully  reconstruct  the  original  input  array  as  long 
as  we  note  the  phase  used  to  obtain  each  particular  resolution  level  from  the  previous 
one.  Such  a  disparity  in  the  sets  of  coefficients  suggests  that  the  DWT  process  of  a  two 
times  decimation  is  not  shift-invariant.  In  other  words,  if  we  decomposed  the  input  array 
with  one  type  of  decomposition  scheme,  our  coefficients  would  be  different  if  the  input 
array  is  shifted  before  the  low  pass  and  band  pass  filtering  processes  occurs.  As  a  note, 
using  the  energy  check  routine,  we  can  see  the  distribution  of  the  energy  of  the 
coefficients  change  as  the  phase  decomposition  scheme  varies.  So  this  shift  variant 
property  can  be  further  exploited  by  varying  the  phase  for  possibly  optimizing  or 
concentrating  the  energy  of  families  of  signals  to  certain  resolution  levels. 

However,  we  desire  linear  shift  invariant  (LSI)  systems  for  the  analysis  of  data. 
For  the  continuous  WT,  this  property  is  generally  true,  however  for  the  DWT  case,  this 
aspect  is  not.  One  way  to  bypass  this  obstacle  is  to  account  for  all  the  phases  during 
each  decomposition.    We  call  this  technique  the  multiple-phase  DWT  (MP/DWT).    A 
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thorough  discussion  for  the  one-dimensional  case  will  first  be  presented,  which  will  then 
be  extended  into  the  two-dimensional  case. 

A.      ONE-DIMENSIONAL  MP/DWT 

Consider  the  case  of  going  from  the  data  input  resolution  m=0  to  m=-l.  The 
coefficients  that  can  be  generated  are  one  of  two  sets,  depending  on  whether  we  selected 
the  phase-0  or  the  phase- 1  case.  We  really  need  only  one  of  the  two  sets  to  decompose 
to  the  next  lower  resolution  level  m=-2  to  maintain  perfect  reconstruction.  Also,  we 
need  to  pick  which  phase  to  decompose  from  m=-l  to  m=-2.  Thus  for  level  m=-2, 
there  are  four  possible  phase  combinations  that  can  result:  00,  10,  01,  and  11.  The 
ordering  of  the  subsequent  phases  is  read  from  left  to  right  to  correspond  to  decomposing 
from  level  m=0  to  m=-2  and  we  call  this  a  phase  vector.  Thus  [iojC.2,5  would  be 
interpreted  as  the  set  of  approximation  coefficients  at  resolution  level  minus  two  and 
decomposed  with  a  phase  vector  of  [10]  with  an  index  of  five.  In  general,  for  any 
resolution  level  m<0,  the  total  number  of  phase  combinations  is  2"m. 

Looking  again  at  going  from  the  m=0  to  m=-l  case  only,  we  can  see  that  by 
sliding  the  filter  coefficients  by  one  instead  of  two  units  to  the  right,  the  resulting 
coefficients  alternate  on  the  phase  selection  starting  with  the  first  coefficient  in  the  phase- 
0  set  followed  by  the  first  coefficient  in  the  phase- 1  set  as  shown  on  the  next  diagram. 
We  determine  the  total  number  of  coefficients  in  level  m=-l  just  like  a  linear  discrete 
convolution,  which  is  the  number  of  data  points  from  level  m=0  plus  the  length  of  the 
filter  mask  minus  one  (  |  {c0}    |  +  N  -  1). 
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Figure  41.   Alternating  of  Coefficients  by  Phase 

Going  from  level  m=-l  to  m=-2  will  require  a  slight  modification  before  shifting 
the  filter  mask  by  one  to  the  right  as  discussed  in  the  previous  paragraph.  Since  two 
sets  of  decomposition  coefficients  are  interleaved,  a  zero  must  be  appended  between  each 
coefficient  in  the  filter  mask  for  properly  decomposing  each  set  of  coefficients.  The  first 
four  values  in  the  resulting  set  correspond  to  the  phases:  00,  10,  01  and  11.  This 
association  repeats  for  each  four  values  in  level  m=-2.  In  fact,  the  spacing  of  each 
coefficient  in  a  desired  phase  is  properly  indexed  by  the  same  amount  as  the  number  of 
zeros  padded  for  display  purposes  in  the  Chapter  III. A. 5. 

This  observation  can  be  generalized  to  any  arbitrary  level  m.  In  this  case,  separate 
the  filter  coefficients  by  2'm-l  zeros  before  the  conducting  single  shift  to  the  right.  By 
noting  the  resolution  level  m,  the  order  of  the  phases  can  be  determined  as  the  binary  bit 
reversal  of  the  binary  values  counting  from  zero  to  2"m-l.  We  show  this  characteristic 
as  a  phase- tree  diagram  for  resolution  levels  m=0  to  m=-3. 
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Resolution  Level 
m  -  0 

m  -  -1 
m  -  -2 


Figure   42 


m  —  -3 


Phase-Tree  Diagram  for  Resolution  Levels  m=0  to 
m=-3 


We  now  define  new  filter  sequences,  h^  and  gm  (m  <0)  for  decomposing  the  values 
Cnm  to  the  next  level  cm.ln.    The  revised  filter  sequences  are  defined  as 


h    =  [  K-N+2)  (OL  K-N+l)  (OL  ...  to)    h(0)  (0)    h(l)  ] 


m  m 


gm  =  [  g(-N+2)  (OL  g(-N+l)  (OL  ...  {OL  *(0)  {0Lm  g(l)  ] 


where  {0}m=2m-l  length  zeros  vector 


(40) 
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Thus,  the  MP/DWT  equations  are  given  as  follows  [Ref.  3:  p.  3][Ref.  4]: 


1 


em-v  -         E         A,<*)V*-i  <41> 

fc=-tf+2-(AM)(2~m-l) 

and 

k=-N*2-(N-l)Q-m-l) 

The  total  number  of  coefficients  also  can  be  generalized  in  terms  of  the  size  of  the 
input  array,  |  do  |  ,  the  size  of  the  h  or  g  mirror  filters,  N,  and  the  current  resolution 
level,  m:  ,.„ 

K\  =  kol  +  (^-D(2w,-i)  <43) 

So  enough  memory  must  be  allocated  if  we  desire  lower  and  lower  resolutions,  since  the 
number  of  possible  phases  increases  by  a  power  of  two  each  step  downward. 

Consider  a  BPSK  signal  shown  in  Figure  43.  Using  a  Haar  wavelet  and  a  phase 
vector  of  [11111111],  a  contour  display  of  the  energy  in  the  approximation  coefficients 
and  of  the  detail  coefficients  results  and  are  shown  next.  Notice  that  much  of  the  180° 
phase  shifts  are  not  readily  detectable. 
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Figure   43.        Sample  BPSK   Signal    [Ref.    4] 
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Figure   44.       Approximation    | c„ 
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loomed      C  o  n  t  o  i_j  r      Display      of     the      Detail  P>Hos«:       11111111 
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Figure   45.       Detail    | d^  | 2 

Now  for  the  MP/DWT  case,  the  180°  phase  shifts  of  the  BPSK  signal  are  now  very 
apparent  in  both  the  approximation  and  detail  coefficients: 


Mu  I  tiphose      Gor-itOLjr-      Disploy      of      the     Approxlmotion 
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Figure    46.        MP/DWT    | c^  | 2  54 
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Figure   47.        MP/DWT    Id^l2 

Twice  the  energy  of  the  coefficients  for  the  current  resolution  level  will  equal  the 
total  energy  on  adjacent  lower  level's  approximation  and  detail  coefficients,  since  there 
are  two  valid  sets  of  coefficients  for  the  reconstruction  to  the  higher  resolution.  So  the 
energy  checking  routine  discussed  in  the  one-dimensional  DWT  is  still  valid  as  long  as 
we  account  for  a  factor  of  one  half  when  going  up  one  resolution  level. 

The  actual  zero  padding  of  the  filter  masks  is  not  the  most  optimum  method  in  the 
determination  of  the  decomposition  coefficients'  proper  phase  order.  A  much  more 
efficient  technique  can  be  determined  by  noting  the  pattern  that  develops  as  Equations 
42  and  43  and  writing  it  out  explicitly.  The  following  diagram  shows  the  developing 
pattern  for  m  =  -2. 
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0  alc(L-2)+aOc(L-4) 

0  alc(L-l  )+aOc(L-3) 

0  0  a0c(L-2) 

0  0  aOc(L-l) 


(N-l) 


m  =  -i 

L  =   |c  0I 
N  =  lh  1=4 


where: 

[aO  al  a2  a3  ...  aN-1] 

[h(-N+2)...h(0)h(l)] 

or 
[g(-N+2)...  g(0)g(l)] 


Figure  48. 


MP/DWT  Pattern  Development  for  m=-2 


Notice  that  the  coefficient  set,  cm_,,  is  a  summation  of  N  vectors  with  |  Cq  |  +  (N- 
l)(2"m-l)  elements.  We  pad  some  of  the  vectors  with  zeros  either  in  the  beginning  or  the 
ending  of  the  approximation  sequence  of  the  current  level  and  multiply  by  an  appropriate 
scalar  value  from  one  of  the  filter  coefficients.  The  following  equations  summarizes  the 
development: 


N-l 


{Vi}  -£*,*.         aP  =  &(-m+p) 


p-0 


p=0 


^  =  figi-N+2+p) 


(44) 


(45) 
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where  *. *  [P-}  W  {24]  (46) 

{Omb}  =  {zeros  coefficient  row  vector  of  size  (N-l-p)*2"m} 
{0^}  =  {zeros  coefficient  row  vector  of  size  p*2"m}. 


This  vectorized  technique  is  faster  than  the  iterative  process  of  padding  the  filter 
coefficients  appropriately.  However,  we  need  more  buffer  space  to  calculate  the  data, 
which  is  dependent  on  the  number  of  resolution  levels  desired. 

The  reconstruction  process  works  in  a  very  similar  fashion  as  in  the  DWT  case. 
All  we  need  is  the  selection  of  the  coefficient  sets,  which  in  turn  determines  the  unique 
path  back  upwards  in  the  phase-tree  diagram,  similar  to  Figure  8.  Once  we  select  the 
desired  phase  at  the  lowest  resolution  level,  the  we  extract  the  corresponding 
approximation  and  detail  coefficients  and  put  them  in  a  packed  format.  Then  the  DWT 
reconstruction  commences  to  get  up  one  level,  where  we  extract  the  next  set  of  detail 
coefficients.   We  repeat  this  process  until  we  reach  resolution  level  m=0. 

B.      TWO-DIMENSIONAL  MP/DWT 

Since  we  define  the  basis  functions  as  separable  in  the  two  orthogonal  directions, 
we  can  use  the  same  arguments  discussed  for  the  one-dimensional  case.  When  we  talk 
about  coefficients  and  masks,  they  are  now  two-dimensional  arrays  as  discussed  in 
Chapter  IV.   Thus,  the  basic  equations  for  the  two-dimensional  MP/DWT  are 
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*w-iW0=2         ££  HhHv(j,k)cm(a+j-l,b+k-\)  m 

;=-^+2-(A//l-l)(2-"-l) 
*=-tfv+2-(Nv-l)(2-m-l) 


I4.-1W»=2        EE  HhGym6m{a+J-\Mk-l)  (48) 

y=-/V2-(ArA-l)(2---l) 
t=-Wv+2-(Wv-l)(2-"-l) 


^m_1(^)  =  2    ££     GhHv(j,k)cm(a+j-i,b+k-D      m 

j=-Nh+2-(Nh-\)(2-m-D 
k—Nv+2-(Nv-l)(2-"-\) 

^m_x{a,b)=2        ££  GhGv(jyk)cm(a+j-hb+k-l)  (50) 

;=-NA+2-(^-l)(2-«-l) 
k=-N+2-(Nv-l)(2-"-\) 


where  the  top  left  corner  of  the  mask  corresponds  to  location  (a,b)  =  (0,0),  and  the 
positive  directions  are  to  the  right  and  downward  [Ref.  3:  pp.  2-3]:. 

Even  the  one-dimensional  vectorized  process  converts  into  a  "matricized"  process. 
Here,  the  resulting  coefficient  array  is  the  matrix  addition  of  NhxNv  total  matrices,  each 
appropriately  padded  with  zeros  on  all  four  sides  and  multiplied  by  one  of  NhxNv  scalar 
values.    The  two-dimensional  equations  result  as  follows: 


M  =   EE"«Z»  ",,  ■  2H,flJL-N*2*p,-N*2*q) 

/>=0     <7=0 


(51) 
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{'<W  =  E  E  *„  z„ 

p=0    ?=0 


b„  =  2HhG(-N+2+p,-N+2+q)  <52) 


pq  h     v 


{2<L,}  •  E  E  <„  z, 

p=0    q=0 


cpq  m  2GhHv(-N+2+p,-N+2+q)  <53> 


p=0     9=0 


f     =  2GhG(-N+2+p,-N+2+q) 


'pq  h     v 


(54) 


This  time  we  define  7^  graphically  in  the  following  figure: 
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Figure  49.  Definition  of  Zn 


zero  coefficient  array 
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The  phase-tree  diagram  can  be  viewed  in  three  dimensions  to  more  appropriately 
show  the  proper  interleaving  of  the  phases.  Since  there  are  four  possible  phases,  there 
will  be  four  more  locations  in  the  indexing  of  the  data  on  the  next  lower  resolution  level. 
The  total  number  of  coefficients  would  be 


[Rows()+(N-l)(2-m-l)iCols()+(Nh-l)(2-m-\)] 


0    V"j| 


(55) 


where  Rows0  and  Cols0  are  the  number  of  row  and  columns  of  the  original  image. 
In  order  to  view  the  phase-tree  diagram  correctly,  we  must  consider  each  resolution  level 
as  a  new  plane  for  the  ordering  of  the  coefficients. 


m  =  0 


m  =  -l 


m  =  -2 


Figure  50. 


Phase-Tree  Diagram 
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The  MP  energy  display  of  the  coefficients  closely  resembles  the  display  for  the 
single  phase  case.  This  energy  display  differs  by  the  fact  that  they  are  the  normalized 
average  energy  values  for  each  level.  As  before,  we  normalize  the  energy  to  the 
resolution  level  m=0.  Consider  as  an  example,  a  MP/DWT  of  the  image  in  Chapter  IV 
using  Haar  wavelets  and  decomposing  down  to  m  =  -3.  The  following  display  shows  the 
MP/DWT  energy  plot: 


wavelet      horlz:      M  a  a  r 
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Figure  51. 


MP/DWT  Energy  Display 
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and  the  energy  check  is  as  follows: 


Level 
m 


Verification  of  the  energy  in  each  resolution  level 

Difference 


Energy  in 
level  m-l= 
c+dl+d2+d3 


Energy  in  C 
level  m 


-2  0.441935911016949 
■1  0.751059322033898 
0      1.000000000000000 


0.441935911016949 
0.751059322033898 
1.000000000000000 


0.000000000000000 
0.000000000000000 
0.000000000000000 


The  MP/DWT  of  the  three-dimensional  energy  display  for  all  the  coefficients  of 
the  image  and  the  corresponding  contour  images  for  the  m=-3  case  are  as  follows: 


wavelet  honz:  Haar 


wavelet  vert:  Haar 


multiphase  case 


Resolution  level  -3 


Figure  52. 


Mesh  Display  of  the  Stored  Energy  for  m  =  -3 
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wavelet  horiz:  Haar 
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100 


50 


50 
multiphase  case 


100 


50 
Resolution  level  -3 


100 


Figure  53.  Contour  Display  of  the  Stored  Energy  for  m=-3 
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VI.    MATLAB  DWT  ROUTINE  DESCRIPTION 

The  Appendix  lists  the  developed  DWT  algorithms.  We  categorize  them  primarily 
into  two  classifications,  the  one  and  two-dimensional  cases.  In  each  situation,  we  further 
define  two  subgroups,  the  single  phase  DWT  and  the  MP/DWT.  Included  with  the 
software  is  an  ASCII  version  of  this  chapter  in  a  README.TXT  file.  This  chapter  first 
describes  the  system  configuration  to  run  the  routines.  Then  we  list  brief  steps  for  each 
of  the  two  general  cases.  These  steps  are  not  totally  inclusive,  but  the  program's 
prompts  will  fill  any  other  gaps  in  the  routine.  Finally,  we  discuss  important 
considerations  for  obtaining  graphical  output  of  the  results  since  this  is  highly  dependent 
on  the  system  hardware. 

A.      SYSTEM  CONFIGURATION 

Minimum  hardware  requirements  for  the  PC  version  are  an  Intel-386 
microprocessor  with  a  math  co-processor,  Microsoft  DOS  or  Digital  Research  DOS  5.0, 
VGA,  at  least  three  megabytes  of  hard  disk  space,  and  four  megabytes  of  RAM,  for  most 
of  the  DWT  applications.  However,  we  recommend  at  least  eight  megabytes  of  RAM 
for  using  arrays  greater  than  10,000  elements,  or  if  m<-5  in  the  two-dimensional 
MP/DWT  case.  Also,  some  DOS  memory  managers  may  conflict  with  the  Pharlap  DOS 
extenders  that  MATLAB  Version  3.5  uses.  For  the  Sun  workstations,  in  order  to  operate 
PROMATLAB  Version  3.5,  we  recommend  the  Open  Windows  or  Sunview  graphical 
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environments  with  at  least  eight  megabytes  of  RAM.  Any  further  memory  can  be 
negotiated  with  the  system  administrator.  Note  that  UNIX-based  operating  systems 
discriminate  characters  between  upper  and  lower  case  letters. 

B.   GENERAL  PROCEDURAL  STEPS 

Initially,  we  must  be  in  the  MATLAB  environment  prior  to  running  any  of  the 
routines.  Either  generate  in  MATLAB  or  load  the  sampled  data  for  resolution  level 
m=0,  and  assign  a  variable  name.  If  and  files  with  the  prefix,  "leg",  and  the  suffix, 
".met"  (all  in  lower-case  letters),  exists  in  the  current  directory,  they  will  be  deleted  once 
the  DWT  routines  commence.  Thus  if  you  desire  to  retain  these  files,  rename  them. 
The  DWT  routines  use  these  files  as  the  output  graphic  files  for  the  routines. 

1.       One-Dimensional  General  Procedures 

a.  If  you  have  no  data,  the  routines  include  three  files  as  data  that  can  be  loaded. 
Type  "load"  followed  by  a  space  and  then  one  of  the  three  names:  "et60.m", 
"et90.m",  or  "bpsk".    The  variable  name  will  be  either  "et60",  "et90",  "bpsk". 

b.  Invoke  "wavld"  in  lower-case  letters. 

c.  Enter  the  variable  name  of  the  input  data. 

d.  Enter  the  sampling  frequency,  <RET>  if  not  known. 

e.  Select  the  desired  h  coefficients. 

f.  Select  "  Y"  if  you  desire  the  MP/DWT  routine  and  skip  to  step  k. 

g.  In  you  did  not  select  the  MP/DWT  routine,  the  single  phase  DWT  commences, 
decomposing  the  coefficients  until  resolution  level  m  =  -  [~log2(  |  cn  |  )"1  with  an 
initial  query  of  if  you  desire  to  see  the  intermediate  resolution  level  displays. 

h.     Choose  either  the  phase-0  or  phase- 1  decomposition  for  each  resolution  level. 
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i.      After  the  final  decomposition,  the  energy  checks  commence. 

j.  Select "  Y"  if  you  desire  to  do  the  reconstruction  process.  If  so,  answer  yes  or  no 
if  you  want  to  display  the  intermediate  levels. 

k.  Lastly,  the  time/scale  display  commences,  and  the  series  of  questions  ask  if  you 
want  to  zoom  in  or  out  and  what  number  of  contour  levels  you  desire. 

1.  For  the  MP/DWT  case,  steps  i  to  k  are  essentially  the  same  except  for  the 
reconstruction,  where  you  must  initially  select  the  desired  phase,  thereby  fixing 
the  phase  path  back  up  to  level  m=0. 

2.       Two-Dimensional  Procedures 

a.  If  you  do  not  have  any  input  data,  invoke  "idata"  for  a  small  selection  of  images. 
The  variable  "im"  (in  lower-case  letters)  will  be  the  input  for  step  c. 

b.  Invoke  "wav2d"  in  lower-case  letters. 

c.  Enter  the  variable  name  of  the  input  data  array. 

d.  Select  the  desired  h  coefficients  in  the  horizontal  direction. 

e.  Select  the  desired  h  coefficients  in  the  vertical  direction. 

f.  Select  "Y"  if  you  desire  the  MP/DWT  routine. 

g.  Enter  the  number  of  resolution  levels  to  go  down.  An  initial  good  selection  is 
four  since  this  will  cover  all  four  possible  phases. 

h.  For  the  single  phase  DWT  case,  select  one  of  four  of  the  decomposition  phases 
for  each  level. 

i.     After  the  final  decomposition,  the  energy  checks  commence. 

j.  For  the  single  phase  DWT  case,  select  "Y"  if  you  desire  the  reconstruction 
routine. 


66 


C.      OBTAINING  GRAPHICAL  OUTPUT 

Except  for  the  initial  time/ frequency  plots,  all  other  displays  will  query  you  if  you 
want  to  store  the  plot.  If  so,  the  plot  is  stored  as  a  metafile  with  the  prefix  of  "leg"  and 
the  suffix  ".met".  A  number  starting  from  zero  is  put  between  the  prefix  and  suffix,  for 
example,  legl4.met  would  correspond  to  the  fifteenth  stored  plot  in  the  DWT  routine. 
Output  of  these  metafiles  are  hardware  and  system  dependent.  Refer  to  the  GPP 
command  in  the  MATLAB  User's  Guide  for  more  detailed  information.  Two 
C-shell  script  files,  called  "metal3"  or  "metal4",  will  generate  a  temporary  Postscript 
file  for  plotting  on  a  Postscript  printer  all  the  metafiles  retained  from  the  DWT  routine. 
The  use  of  these  files  is  currently  only  guaranteed  to  work  at  the  US  Naval  Postgraduate 
School  Electrical  Engineering  Department's  Sun  workstations,  and  we  list  them  here  as 
a  guide  in  generating  other  C-shell  script  files  for  a  particular  system.  For  the  DOS 
version,  the  following  command  will  generate  all  the  metafiles  onto  a  HP  Laserjet  III 
printer  as  long  as  the  we  invoke  the  following  DOS  command  in  the  same  directory  as 
the  metafiles:  for  %f  in  (leg*. met)  do  gpp386  %f  /djetl50  /ol  /fprn  .  If  you  put  it  in 
a  DOS  batch  file,  echo  the  percent  sign  (  %%  instead  of  %  ).  Please  note  that  the 
computer  may  take  some  time  generating  these  plots.  Be  also  aware  that  depending  on 
the  size  of  the  plot,  the  graphic  files  may  exceed  the  printer  buffer  memory. 
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VII.    CONCLUSIONS 

In  this  thesis,  we  developed  discrete  wavelet  transform  algorithms  for  both  the  one 
and  two-dimensional  cases.  The  iterative  mechanism  of  decomposing  the  data 
coefficients  by  a  convolution-and-subsample  process  is  achievable.  These  realizations 
require  that  the  scaling  and  wavelet  functions  form  orthonormal  sets  with  compact 
support.  Also  for  the  two-dimensional  case,  the  basis  functions  were  separable  in  the 
two  orthogonal  directions. 

We  found  that  the  output  of  the  decomposition  coefficients  can  be  causal  provided 
that  the  definitions  of  the  indices  for  both  filter  coefficients,  h  and  g,  must  be  selected 
to  accommodate  the  causal  condition.  In  fact,  the  definition  of  the  g  filter  coefficients 
must  be  adjusted  from  the  definitions  given  by  Mallat  and  Daubechies.  We  found  for 
the  one-dimensional  case,  the  decomposition  is  simply  a  dot  product  of  the  filter 
coefficients  with  a  set  of  values  windowed  from  the  higher  resolution  level  approximation 
coefficients.  This  window  subsequently  shifts  two  units  to  the  right  for  the  next 
decomposed  coefficient.  We  found  that  the  reconstruction  process  uses  the  same 
decomposition  algorithm  if  we  separate  each  of  the  individual  input  values  by  a  zero,  the 
filter  coefficients  are  time-reversed,  and  the  shift  is  now  one  unit  to  the  right  at  a  time. 

The  two-dimensional  case  was  a  natural  extension  of  the  one-dimensional  process. 
We  found  that  the  decomposition  can  be  realized  as  a  mask  element  multiplication  and 
summation,  with  the  mask  shifting  two  units  to  the  right  until  the  row  completion  and 
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then  two  units  down  to  process  the  next  row  of  coefficients.  The  shifting  of  the  mask 
is  similar  to  a  raster  scan.  The  reconstruction  follows  the  same  process  as  in  the 
decomposition  but  each  of  the  input  array  coefficients  are  separated  by  a  row  and  a 
column  of  zeros,  the  filter  masks  are  spatially  reversed  in  both  direction,  and  the  raster 
shift  is  one  unit  instead  of  two. 

We  found  that  there  were  two  possible  phases  in  decomposing  the  data  in  the  one- 
dimensional  case,  and  four  possible  phases  for  the  two-dimensional  case.  This  led  to  the 
conclusion  the  current  decomposition  algorithm  was  very  shift  variant,  contrary  to  the 
continuous  version  to  the  wavelet  transform,  which  is  shift  invariant.  We  then 
determined  that  the  discrete  case  can  possibly  approximate  the  continuous  version  if  we 
account  for  all  the  phases  during  the  decomposition.  We  then  developed  the  multiple- 
phase  discrete  wavelet  transform  for  both  one  and  two  dimensions.  The  coefficient  sets 
for  a  particular  phase  can  be  interleaved  with  other  sets  of  a  different  phase.  We  found 
that  the  multiple-phase  was  more  capable  in  detecting  discontinuous  properties  of  data 
than  in  the  regular  discrete  wavelet  transform,  such  as  a  binary  phase  shifting  function. 
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APPENDIX  DWT  MATLAB  FILE  LISTINGS 


A.      ONE-DIMENSION AL  ROUTINES 


%  15  SEP  92 

%  WAV1D.M        The  One  Dimensional  Discrete  Wavelet  Transform 

% 

% 

%  By:       LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

% 

% 

%  This  routine  does  a  one-dimensional dicrete  wavelet  decomposition 

%  of  a  one-dimensional  data  array.    The  algorithm  allows  for  the 

%  selection  of  the  desired  phase  for  the  respective  decompositon 

%  level.    The  following  routines  are  necessary  for  the  proper 

%  operation  of  this  algorithm: 

% 

%     chkputer.m  --  algorithm  checks  for  compatible  hardware 

%     daubdata.m  -  determines  the  h  coefficients  for  a  daubechies 

%  compactly  supported  orthogonal  wavelet 

%     mp.m  -  MultiPhase  decomposition 

%     phsl.m        --  selects  the  phase- 1  decomposition 

%     phsO.m       -  selects  the  phase-0  decomposition 

%     dpscoef.m  —  display  the  coefficients  for  each  level 

%     enrgld.m    -  determines  the  energy  check  for  each  resolution 

%     reconld.m  —  reconstruction  algorithm  for  verification 

%     plane. m      —  displays  the  coefficients  in  the  time-scale  domain 

%     strplt.m     -  stores  the  plots  in  meta-files 

% 

clc 

chkputer    %  Checking  for  the  basic  hardware  necessary 


% 

%  cO  corresponds  to  the  resolution  level  0,  the  original  data 

%      array. 

clc 

c0  =  input(' Enter  the  name  of  the  data  in  vector  row  format:  '); 

%  We  will  check  if  the  data  is  in  row  format,  if  in  column  format, 
%  take  the  transpose  of  the  input. 

[il  jl]=size(c0); 

if         jl  =  =l,cO=cO'; 
elseif     il>l&jl>l, 

disp('Bad  Input  Data,  Restart  the  Program') 
return 
end 
clear  il  jl 
% 
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fsamp  =  input('Enterlhe  sampling  frequency  (Hz):    '); 


% 

%  The  following  code  determines  the  desired  wavelet 

for  i  =  l:8,disp(['  ']),end 

ql  =  ['Enter  the  number  for  the  desired  choice:  ' 


(1)  Haar  h  coefficients 

(2)  Daubechies  h  coefficients 

(3)  Own  set  of  h  coefficients 


disp(ql) 

qqq  =  input('  '); 

if         qqq==3, 

h  =  input(' Enter  your  own  set  of  h  coefficients'); 
%  Checking  if  the  data  are  compact  support 

a=sum(h);b=h*h'; 

toll  =  le-12; 

while  abs(l-a)  >loll  &  abs(.5-b)  >  toll 
disp('Yourh  coefficients  are  NOT  compactly  supported!') 
h  =  input('Re-enterlhe  h  coefficients  or  ctrl-z  to  stop'); 

end 

wwavlet  =  'OTHER'; 
elseif   qqq==2, 

qq  =  input('Howmany  taps  (2-10)?  '); 

h=daubdata(qq); 

wwavelel  =  ['Daubechies  ',num2str(qq),'-Up']; 
else      h  =  [.5  .5];wwavelet='HAAR'; 
end 

h  =  sqrt(2)*h;      %  The  factor  of  sqrt<2)  normalizes  the  energy 

% 

%  Some  definitions  of  variables  used  throughout  the  entire 
%         algorithm: 

% 

%  wwavelet  —  name  of  the  selected  wavelet 

%  LL  --the  length  of  the  input  array 

%  Nh  -  the  number  of  h  and  thus  g  coefficients 

%  lowest     —  the  number  of  resolution  levels  below  the  input 

%  pltcnt     —  plot  counter  for  storing  desired  plots 

%  h  —  the  "h"  coefficients 

%  g  —  the  'g*  coefficients 

%  lowest     -  the  lowest  resolution  level 

%  phsvet    —  the  phase  vector  for  recording  the  appropriate  phase 

%  zro         --  zero  vector  to  record  the  appropriate  zero  padding 

%  for  each  resolution  level. 


%  We  will  intialize  some  of  the  constants: 

LL=length(cO); 

Nh  =  length(h); 

lowest  =  -ceil(log(LL)/log(2)); 

lvl  =  0; 

pltcnt  =  0; 

% 
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%    Generating  the  g  coefficient* 

g  =  fliplr(h); 
fori=2:2:Nh 

g(l,i)  =  -g(l.i); 
end 
% 


% 

%  The  following  set  of  lines  branches  to  the  multiphase 

%  decomposition  algorithm  if  desired.    Program  ends  after 

%  the  multiphase  execution 

clc 

q  =  input('Do you  desire  the  muliphase  decomposition  (Y/N)?  '.'»')', 

ifq=  =  'Y'  |  q=  =  Y 

mp 

return 
end 
% 


% 

%  The  Decomposition  Routine 

%  Determine  if  we  want  to  ignore  the  display  for  each  resolution 

%        level's  coefficients. 

pltl  =inputCBypas8the  display  of  the  Coefficients  (Y/N)?  \V); 

clc 

phs  =  2;    %  This  just  sets  "phs"  initially  to  be  NOT  equal  to  1  or  0 

phsvct  =  []; 


%  

%  The  following  code  sets  the  number  of  variables  necessary 
%  to  record  the  coefficients,    ie.    lowest  =  -2.  we  have 
%  cO,cl,c2,dl,d2 

1.11  =1.1  ; 

for  lvl  =  -l:-l:lowest 
while  1  =  =  1 
phs  =  input(' Enter  the  desired  Phase  [0/1]  for  this  resolution  level:  '); 
ifphs=  =  l  | phs  =  =0, break, end 
end 

phsvct  =  [phs  phsvct]; 
eval([  'c '  ,num2str(ab8(lvl)) , '  =  phs '  ,num2stx(phs) ,  '(c' , . . . 

num2str(ab8(rvl)-l),\h);']) 
eval([  'd '  ,num2str(ab8(lvl)) , '  =  phs'  ,num2str(phs) ,  '(c' , . . . 
num2str(abs(lvl)-l),',g);']) 

%  This  if  statement  runs  the  display  of  the  coefficients  if  the  flag 
%         pltl  is  set  to  NO 
ifpltl  =  =  'N'  |  pltl  =  =  'n',  dspcoef.  end 
phs  =  2; 

eval(['LLL=rnax(LLL,length(c',num2str(-lvl),')*2*(-lvl));']) 
zro  =  Izro2*(-lvl)-lJ; 
end 


%  We  now  do  an  energy  check  of  the  coefficients 
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enrgld 
clg 


% 

%  We  now  run  the  reconstruction  option 
reconld 


% 

%  We  now  display  the  coefficients  in  the  time-scale  domain 

plane 

% 

%  END   OF  WAV1D.M 


% 

%  Phase  Zero  decomposition  for  the  1-D  case 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

% 

%  phsO.m  is  a  function  m-file  the  does  a  phase  0  decomposition  of 

%  of  the  data  array  using  either  the  h  or  g  coefficients.   The  main 

%  program  that  calls  this  routine  is  "wavld.m" 

% 

%  The  input  arguments  for  this  routine: 


% 

data 

-    input  data 

% 

h 

-   the  h  or  g  coefficients 

% 

L 

-   the  length  of  the  data  vector 

% 

H 

-   the  length  of  the  h  or  g  coefficient  vector 

function  x=phsO(data.h) 

L=length(data); 
H  =  length(h); 


%  The  following  4  lines  check  if  the  data  array  is  even,  for  the 
%        phase  0  case,  we  must  pad  it  with  one  zero  if  true 
if  rem(L/2,floor(L/2))  =  =  0 
data  =  [data  0]; 
L=L+1; 
end 

data  =  [zeros(l,H-l)datazeros(l,H-2)]; 
L=L+2*H-3, 
a  =  0; 
forii  =  l:2:L-(H-l) 

a  =  a+l; 

xa.a^dataGKii  +  H-iriT; 
end 
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% 

%  Phase  One  decomposition  for  the  1-D  case 

% 

%  By:      LT   J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece. nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

% 

%  phsl  .m  is  a  function  m-file  the  does  a  phase  1  decomposition  of 

%  of  the  data  array  using  either  the  h  or  g  coefficients.    The  main 

%  program  that  calls  this  routine  is  'wavld.m* 

% 

%  The  input  arguments  for  this  routine: 

% 

%  data     --    input  data 

%         h         —   the  h  or  g  coefficients 

%         L         --   the  length  of  the  data  vector 

%         H  the  length  of  the  h  or  g  coefficient  vector 

% 

function  x  =  phsl  (data, h) 

L  =  length(data); 
H=length(h); 

%  The  following  four  lines  checks  if  the  data  vector  is  odd,  if  it 
%        true  for  the  phase  1  case,  we  must  pad  it  with  one  zero 
if  rem(L/2,floor(L/2))  ~  =  0 

data  =  [data  0]; 

L=L+1; 
end 

data  =  [zeros(  1  ,H-2)  data  zeros(  1  ,H-2)] ; 
L=L+2*H-4; 

a  =  0; 

forii  =  l:2:L-(H-l) 

a  =  a+l; 

x(l,a)=dato(l,ii:ii  +  H-l),h'; 
end 


% 

%  dspcoef.m      Display  the  coefficients 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School.  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.np8.navy.mil 

%  Phone:     (408)  646-3044/2772 

% 

%  dspcoef.m  displays  the  coefficients  for  each  particular 

%  resolution  level  with  the  appropriate  scaling  factor.    There  is 

%  also  a  special  routine  embedded  in  this  algorithm  for  the 

%  proper  bar  display  of  the  HAAR  case.  The  main  program  is  "wavld.m" 

% 

%  The  variables  are  defined  as  follows: 

% 

%  z  -  the  number  of  zeros  to  be  padded  btwn  coeffs 

%  lvl         --  variable  from  "w.m"  for  resolution  level 
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%  fsamp  -  sampling  frequency  of  the  data 
%  wwavelet  --  string  for  the  selected  wavelet 
* 

clg 

hold  off 
axis( 'normal') 

%  For  the  HAAR  case  do  the  following: 
if  wwavelet(  1 )  ==  'H'   |  wwaveletf.1)  ==  'h' 
n=2*(-lvl); 

eval(['ctemp  =  2A(lvl/2)*c',num2su-(-lvl),y]) 
eval(['dtemp  =  2*(lvl/2)^.num2slr(-lvl),';']) 
ymax  =  max(max(ctemp),max(dtemp)); 
if  ymax<0,ymax  =  0;end 
ymin  =  min(min(clemp),min(dtemp)); 
if  ymin>0,ymin  =  0;end 
L  =  length(ctemp) ; 
if  L=  =  l 

axis([0  (n*1.5-l)/fsamp  1.2*ymin  1.2*ymax]); 

subplot(211) 

plot([0  0  n*1.5-l  n*1.5-l]/fsamp,[0ctemp  ctemp  0]) 

tiUe([ 'Approximation at  Resolution  level  ',num2slr(lvl)]) 

xJabel('Haar  wavelet') 

axis([0  (n*1.5-l)/fsamp  1.2*ymin  1.2*ymax]); 

subplot(212) 

plot([0  0  n*1.5-l  n*1.5-l]/fsamp,[0dtemp  dlemp  0]) 

tiUe<['Detail  Phase- ',num2slr(phs)]) 

pause 

itiplt 

else 
x=[0:n:n»H]  +  .5»n; 

x  =  x/fsamp; 

axis([0  x(length(x))  1.2*ymin  1.2*ymax]); 
subplot(21  l),bar(x, ctemp) 

u'Ue([ 'Approximation at  Resolution  level  '.num2slr(lvl)]) 
xlabel('Haar  Wavelet') 
axis([0  x(length(x))  1.2*ymin  1.2*ymax]); 
subploi(212),bartx.ditmp) 

UUe<[' Detail  Phase-  ',num2str(phs)]) 

pause 
strplt 
end 

%  For  NON-HAAR  cases,  do  the  following: 
else 

z  =  2*(-lvl)-l; 

eval([  'ctemp  =  2*(lvl/2)*ptzero(c '  ,num2slr(-lvl) , '  ,z) ; ']) 

eval(['dtemp  =  2*(lvl/2)*ptzero(d',num2str(-lvl),',z);'J) 

ymax  =  max(max(ctemp)  ,max(dlemp)) ; 

if  ymax<0,ymax  =  0;end 

ymin  =  min(min(  ctemp)  .min(dtemp)); 

if  ymin  >0, ymin  =  0;end 

%  We  can  display  the  data  with  a  window  size  of  the  original 
%       input,  or  show  the  spillover  values 
q  =  ['Enter  the  type  of  display:' 


•    1  - 

Windowed 

Display 

'     2- 

Non-  Windowed 

clc 

disp(q) 

disp(f 

']') 

q  =  input('  '); 
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if  q=  =2 

x  =  [0:length(ctemp)-l]  +  .5; 

x  =  x/fsamp; 

axis([0  x(lenglh(x))  ymin  ymax]); 

subplot(21  l),bar(x,clcmp) 

tille([ 'Approximation at  Resolution  level  ',num28tr(lvl)]) 

xlabel([wwavelet,'  Wavelet']) 

axis([0  x(lenglh(x))  ymin  ymax]); 

subplot(2 1 2)  ,bar(x,dl£mp) 

title([' Detail  Phase- ',num2str(phs)]) 

xlabel(['Non-Windoweddisplay    -  max  windowed  value:  ',... 
num2str(length(x)- 1 )]) 

pause 

strplt 
else 

x  =  [0:LL-l]  +  .5; 

x  =  x/fsamp; 

axis([0  x(length(x))  ymin  ymax]); 

subplot<2 1 1)  ,bar(x,ctemp(  1,1:  LL)) 

titled 'Approximation  at  Resolution  level  ',num2str(lvl)]) 

xlabel([wwavelet,'  Wavelet']) 

axis([0  x(length(x))  ymin  ymax]); 

subplot^  1 2).bar(x,dtemp(l ,  1  :LL)) 

title(['Detail  Phase- ',num2str(phs)]) 

xlabel(['Windowed display  —  max  windowed  value:  ',num2str(LL-l))) 

pause 

strplt 
end 
end 


%  enxgld.m        Determining  the  Energy  in  each  resolution  level 

%  One  Dimensional  Case 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

% • • •••• • 

%  For  the  single  phase  case,  the  energy  is  normalized  to  resolution 

%  level  0.    AU  the  energies  of  the  'cV  and  the  "d's"  of  a 

%  particular  level  should  add  up  to  the  energy  of  the  "c's"  of  the 

%  next  higher  level. 

% 

%        enrgc  -  energy  array  of  the  'c's* 

%        enrgd  —  energy  array  of  the  "d's" 

%        norm   —  energy  of  the  data  "cO" 

%        esum   -  energy  sum  array  (enrgc  +  enrgd) 

% 

%  The  calling  routine  is  "wavld.m"  . 

%  strplt.  m  is  a  necessary  m-file 


enrgc =zeros(  1 ,  -lowest); 
enrgd  =  zeros(  1 ,  -lowest) ; 
a  =  0; 

norm  =  sum(c0  .*2); 
for  ii  =  lowest:  1:-1 

a  =  a+l; 

eval([  'enrgc(a)  =sum(c"  ,num28lr(-ii) , ' .  "7) ; ']) 
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eval([  'enrgd(a)  =  8um(d '  ,num2slr(-ii), ' .  *2); ']) 
end 

enrgc  =  enrgc/norm; 
enrgc  =  [enrgc  1]; 
enrgd  =  enrgd/norm, 
clg 

axis([lowest-2  1  0  1.2]); 
subplot(2 1 1 )  ,baK[lowest:  0] ,  [enrgc]) 
tille(['Normalized  energy  of  the  "c"  coefficients    Phase; 

setstr(flipu-(phsvct+48))]) 
xlabel('Resolution  level') 
axis([lowest-2  1  0  1.2]); 
subplot<2 1 2)  ,bar([  lowest:  - 1  ] ,  [enrgd]) 
tille(' Normalized  energy  of  the  "d"  coefficients') 
xlabclf.' Resolution  level') 
pause 
strplt 


%    Energy  check  of  the  previous  plot 

esum  =  enrgc  +  [enrgd  0] ; 

clc 

dd  =  ['         Verification  of  the  energy  in  each  resolution  level 


'Level        Energy  in  Energy  in  C  Difference' 

level  ml  level  m 

c  +  d 

]; 

disp(dd) 

j  =  -lowest; 

for  i=  1: -lowest 

a  =  e8um(i);b  =  enrgc(l,i+  l);c  =  l-j;d  =  abs(a-b), 
fprintf('%2.0f       %17.15f         %17.15f      '.c.a.b) 
fprintf('%17.15f\n\d) 

j=j-i; 

end 
pause 
clc 
axis; 


% 

%  reconld.m    One  Dimensional  Single  Phase  Reconsruclion 

% 

%  By:       LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

% 

% 

%  This  routine  conducts  a  reconstruction  of  the  the  data  with  the 

%  proper  phase  selected  between  resolution  levels.    This  is  a  sub- 

%  routine  for  the  main  wavelet  program  "wavldm". 

% 

%  The  following  variables  are  farther  defined: 

% 

%        hh  &  gg  -  reconstruction  coefficients 

%        lowest  -  defined  in  "wavldm" 

%        cwork  —  reconstructed  c'values  at  a  particular  lvl 

%        phsvct  -  phase  vector  defined  in  "wavldm". 
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% 

%  The  following  m-files  are  necessary  for  proper  operation: 

% 

%  wavld.m     --  main  program 

%  rphsl.m     --  reconstruction  phase- 1  function  m-  file 

%  rphsO.m     --  reconstruction  phase-0  function  m-file 

%  strplt. m   -- storage  of  the  displays 


clc; 

disp(['  •]') 

q  =  input('Doyou  desire  to  do  the  Reconstruction  Comparison  (Y/N)?',V); 

clc 

ifq  ==  -V   |  q  ==  y, 
%  generation  of  the  reconstruction  coefficients 
hh=fliplr(h);gg  =  fiiplr(g); 
ii  =  1 ; 

eval([  'cwork  =  c'  ,num2str(-lowest) ,';']) 
for  lvl  =  lowest:  1:-1 

eval(['dwork=d',num2str(-lvl),';']) 

eval([ 'cwork  =  rphs',num28tr(phsvct(ij)),'(cwork,d  work, hh.gg);']) 

eval(['a  =  length(c',num28tr(-lvl-l), ');']) 

cwork  =  cwork(  1 :  a) ; 

ii=ii+l; 

phswrk=phsvct(l:length(phsvct)-l); 

%     Comparison  of  the  Original  and  the  Reconstructed  data 

clg 

subplot(211) 

eval([•bar(c•,num28tr(-lvl-l),')•]) 

utle([ 'Original  Data        Phase:  ',setstr((phswrk+48))]) 

xlabel([ 'Wavelet:  \wwavelet]) 
subplot(212) 

bartcwork); 

tit]e( 'Reconstructed  Data') 

xlabel(['Resolulion  level  ',num2str(lvl+l)]) 
pause 
strplt 
clg 

a  =  eval(['c',num2su-(-lvl-l),'-cwork']); 
bar(abs(a)); 

title( 'Absolute  Error  in  the  Reconstruction') 
xlabel(['Re«oluuon  Level  ',num2str(lvl+ 1)]) 
pause 
strplt 
clg 

phswrk  =  phswrk(  1 :  length(phswrk)- 1 ) ; 
end 
end 


function  out  =  rphsO(cwork,dwork,hh,gg) 

%***••***•••••»••*••*•*••••*•••»*•••**•************••»*»**»••***•*** 

%  rphsO.m      Reconstruction  Phase  Zero  routine 

% 

%  By:  LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 
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%  This  routine  does  a  reconstruction  of  a  particular  level  with  a 

%  phase  0  selection.    The  calling  program  is  "wavld.m".    Input 

%  arguments  are: 

% 

%  cwork   —    "c"  coefficients  of  a  particular  level 

%  dwork   -    "d"  coefficients  of  a  particular  level 

%  hh        —    reconstruction coeffs  for  the  "c's" 

%  gg        -    reconslruciton  coeffs  for  the  "d's" 

% 

%  "out"  is  the  output  argument!  that  equals  the  c  coeffs  of  the 

%    next  higher  resolution  level. 

% 

L  =  length(c  work) ; 
xl=D; 
x2  =  D; 
forj  =  l:L 

xl=[xlcwork(l,j)0]; 

x2  =  [x2dwork(l,j)0]; 
end 

xl=[xl  zeros(l.length(hh)-2)]; 
x2  =  [x2  zeros(l,lenglh(gg)-2)]; 
a  =  0; 
for  j  =  1  :length(xl)-(length(hh)- 1) 

a  =  a  +  l; 

out(l,a)=xl(l,j:j  +  length(hh)-l)*hh'  +  x2(l,j:j  +  length(gg)-l),gg' 
end 


function  out  =  rphsl (cwork, dwork, hh.gg) 

% • ......................... 

%  rphsl.m      Reconstruction  Phase  One  routine 

% 

%  By:  LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

Or******************************************************************* 

% 

%  This  routine  does  a  reconstruction  of  a  particular  level  with  a 

%  phase  1  selection.    The  calling  program  is  "wavld.m".    Input 

%  arguments  are: 

% 

%  cwork   -    "c"  coefficients  of  a  particular  level 

%  dwork    --    "d"  coefficients  of  a  particular  level 

%  hh       —    reconstruction  coeffs  for  the  "c's" 

%  gg        --    reconslrucilon  coeffs  for  the  "d's" 

% 

%  "out"  is  the  output  argumemt  that  equals  the  c  coeffs  of  the 

%    next  higher  resolution  level. 

% 

L  =  length(cwork) ; 
xl=[|; 
x2  =  []; 
forj  =  l:L 

xl=[xl  cwork(l,j)0]; 

x2  =  [x2dwork(lj)0]; 
end 

xl  =[0  xl  zeros(l,length(hh)-2)]; 
x2  =  [0  x2  zeros(l,length(gg)-2)]; 
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«  =  0; 

for  j  =  1 :  length(x  1  )-(length(hh)- 1 ) 

a  =  a  +  l; 

oul(l,a)  =  xl(l,j:j  +  length(hh)-l)*hh,  +  x2(l,j:j  +  lenglh(gg)-l)*gg-, 
end 


«»•*••»•*»»•»••*•••*♦•*•»*»•*•*••••********•***••********»•**»*••••* 

%  plane. m     develop  the  phase  plane  diagram 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.np8.navy.mil 

%  lam@ece.np8.navy.mil 

%  Phone:     (408)  646-3044/2772 

or******************************************************************* 

%  Phase  Plane  Diagram  for  the  One  Dimensional  Single  Phase  case 

%  This  routine  is  called  by  the  parent  routine  "wavld.m".   All  the 

%  variables  defined  in  'wavld.m'  are  used  here.  New  variables: 

% 

%       c  (1-lowest)  X  LL  "c"  coefficient  matrix 

%      d  (-lowest)  X  LL  "d"  coefficient  matrix 

%       N.N1  —  variables  for  the  number  of  contour  levels 

% 

%  The  main  subroutine  used  is  ptzero.m  funtion  m-file 

% 

clc 
clg 

%    Set  up  the  c  and  d  matrices  from  the  "c*"  and  "d*"  variables 

c  =  zeros(-iowe8t+  1  ,LL); 
d  =  zeros(-lowest,  LL) ; 
c(l,l:length(cO))  =  cO; 
clear  cwork  dwork 
for  lvl=  l:-lowest 

cwork  =  eval([  'ptzero(c '  ,num2str(lvl) , '  ,zro(lvl)) ']); 

dwork  =  eval([  'ptzero(d '  ,num2str(lvl), '  ,zro(lvI)) ']); 

c(lvl  +  1,1:  length(cwork))  =  cwork; 

d(lvl,  1 :  length(dwork))  =  dwork; 
end 

c  =  flipud(c.*2); 
d  =  flipud(d.*2); 
[i  j]=size(c); 

x  =  0:j-l;yc  =  0:i-l;yd  =  l:i-l; 
mesh(c); 
title([' Energy  Display  of  the  Approximation    Phase:  ',... 

setstr(fliplr(phsvct+48))]); 
pause 
N  =  10; 

contour(c,10,x,yc); 

title<[' Contour  Display  of  the  Approximation    Phase:  ',... 
8etstr(fliplr(phsvct+48))]); 

ylabel('-n  resolution  level'); 

xlabel('10  contour  levels') 
pause 
clg 

ql  =input( 'Do you  desire  to  zoom  in  on  the  display  (Y/N)?',V); 
whileql  =  =  'Y'  |  ql   ==  'y', 

q2  =  input('Doyou  want  to  see  the  plots  again  (Y/N)?',V); 

ifq2  =  =  'Y'  |  q2  ==  'y', 

mesh(c);title<'Previous  Energy  Display  of  the  Approximation'); 
pause 
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contour(c,N,x,yc);uUe<'Previou8Contour  Display  of  the  Approximation'); 

ylabel('-n  resolution  level'); 

xlabel(['Numberof  counlours:  ',num2Blr(N)]) 

pause 
end 
clc 

disp(['     ']') 

disp(['x  range:  0  to  ',num2str(j-l)]) 
disp(['y  range:  0  to  -',num2str(i-l)]);di8p([       ]') 
xl  =input( 'Enter  the  minimum  x-value:')+  1; 
x2  =  input(' Enter  the  maxumum  x-value:')+  1; 

yl  =input('Enlerthe  minimum  'magnitude'  y-value  (least  negative): ')  +  l; 
y2  =  input{'Enterthe  maxumum  'magnitude'  y-value  (most  negative): ")+l; 
yl  =abs(yl);y2  =  abs(y2); 
y  =  l:l:i;y  =  fuplr(y); 
yll=y(y2); 
y22=y(yl); 
clg 

me8h(c(yll:y22,xl:x2)) 
title(['Zoomed  Energy  Display  of  the  Approximation     Phase:  '.... 

setstr(fliplr(phsvct+48))]); 
pause 
strplt 
clg 

Nl  =  10; 
while  Nl<  999 
N  =  N1; 

contour(c(yll:y22,xl:x2),Nl,xl:x2,yl-l:y2-l) 
title([' Zoomed  Contour  Display  of  the  Approximation     Phase:  ',... 

aetstr<fliplr(phsvct + 4  8))]) ; 
ylabel('-n  resolution  level') 
xlabel(['Number  of  contours:  ',num2str(Nl)]) 
pause 
strplt 

Nl  =input( 'Enter  the  number  of  contour  levels  (999  to  continue)  '); 
end 

ql  =  input('Zoomin  or  out  Further  (Y/N)?  \V); 
end 
clc 

%    Now  the  Detail  Signal 

mesh(d); 

title(['Energy  Display  of  the  Detail      Phase:  \setstr(fliplr<phsvct+48))]); 

pause 

clg 

contour(d,10,x,yd); 

title(['Contour  Display  of  the  Detail     Phase:', setstr(fliplr(phsvct+48))]); 

ylabel('-n  resolution  level'); 

xlabel('10  contour  levels') 

pause 

clg 

clc 

ql  =input('Doyou  desire  to  zoom  in  on  the  display  (Y/N)?',V); 

while  ql  =  ='Y'  |  ql  =  =  'y', 

q2  =  input('Doyou  want  to  see  the  plots  again  (Y/N)?V>'); 
ifq2  =  =  'Y'  |  q2  =  =  'y'. 
mesh(d);title('Previous  Energy  Display  of  the  Detail'); 
pause 

contoured, N,x,yd);title('PreviousContour  Display  of  the  Detail'); 
ylabel('-n  resolution  level'); 
xlabel(['Number  of  contours:  ',num2str(N)]) 
pause 
end 
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disp(['x  range:  0  to  ',num2str(j-l)]) 

disp(['y  range:  -1  to  -',num2str(i-l)]);disp([       ]') 

xl  =  input('Enterthe  minimum  x-value:')+  1; 

x2  =  input( 'Enter  the  maxumum  x-value:')  +  l; 

yl  =input('Enterthe  minimum  "magnitude"  y-value  (least  negative):'); 

if  yl  ==  0,  yl=l;  end; 

y2  =  input(' Enter  the  maxumum  "magnitude"  y-value  (most  negative):'); 

yl  =abs(yl);y2  =  abs(y2); 

y  =  l:l:i-l;y  =  fliplr(y); 

yll=y(y2); 

y22=y(yl); 

clg 

mesh(d(yll:y22,xl:x2)) 

title([ 'Zoomed  Energy  Display  of  the  Detail     Phase:  ',... 

setstr(fliplr(ph8vcl + 48))]) ; 
pause 
strplt 
clg 
while  N<  999 

conlour<d(yll:y22,xl:x2),N,xl:x2,yl:y2) 

tille(['Zoomed  Contour  Display  of  the  Detail    Phase:  ',... 
setstr(fliplr(phsvct+48))]); 

ylabel('-n  resolution  level') 

xlabel([ 'Number  of  contours:  ',num2str(N)]) 

pause 

strplt 

N  =  input(' Enter  the  number  of  contour  levels  (999  to  end):  '); 
end 

ql  =input('Zoomin  or  out  further  (Y/N)?  W); 
end 


%    Multiphase  routine  for  The  Discrete  Wavelet  Decomposition  Routine 

% 

%  By:       LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

% 

% 

%    This  routine  does  a  multiphase  discrete  wavelet  decomposition  of 

%    a  one-dimensional  data  vector.    This  routine  requires  the  following 

%    M- files: 

%  wavld.m 

%  strplt.  m 

%  zOpad.m 

% 

maxSamp  =  (Nh-l)»(2'(-lowe8t)-l)  +  LL; 
coeffc  =  zeros(- lowest  +  1  .maxSamp); 
coeffc(-lowest  +1.1:  LL)  =  cO; 
coeffd  =zeros(-lowest.maxSamp); 
numcoef =zeros(-lowest+ 1,1); 

% 

% 

%    ***   The  Main  Multi- Phase  Decomposition  Algorithm   *** 

numcoef(-  lowest  +  1 )  =  LL; 

>  =  1; 

for  row  =  -lowe«t:-l:l 

numcoef(row)  =  (Nh-l)*(2*i-l)  +  LL; 

pi  =zOpad(coeffc(row+l,l:numcoef(row+l)),h,i); 
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coeffc(row.  1 :  lenglh(pl ))  =  pi ; 

p2  =  zOpad(coeffc(row+  l,l:numcoef(row+  l)).g,i); 
coeffd(row,  1  :length(p2))  =  p2; 
clear  pi  p2 
i  =  i+l; 
end 

% 

%        Display  of  Coefficients  at  each  resolution  level 
clc 

q  =  input('Do  you  desire  to  see  values  at  diff  resolution  levels  (Y/N)?  \V); 
ifq  ==  'V   |  q==  y 
plotr.^LL-lLcO,'*') 

xlabel(['n  sampling  points  n*',num2str(l/fsamp). 'seconds']); 

ylabel( 'Amplitude') 

title( 'Original  Data  (resolution  level  0)    Multiphase  case'); 
pause 
strplt 
clg 

index  =  [-lowest:- 1 : 1] ; 
xl=0;x2  =  -lowest+l; 
clc 

disp(['Levels  to  pick  are  from  -1  to  '.num2str<lowe8t)]) 
while  xl<l   |  x2> -lowest 
xl  =input(' Enter  first  resolution  level  of  the  desired  range:  '); 
x2  =  input(' Enter  second  resolution  level  of  the  desired  range:  '); 
xl  =abs(xl);x2=abs(x2); 
if  xl  >x2,  temp  =  x2;  x2  =  xl;  xl  =lemp;  end 
end 
while  xl  <  =x2 

pts  =  numcoef(index(xl)); 
Tc  =  coeffc(index(x  1 ) ,  1 :  pts) ; 
Td  =  coeffd(index(xl).  1  :pts); 

minpts  =  min(min(Tc)  ,min(Td)) ;  if  minpts  >  0 ,  minpls  =  0 ;  end 
maxpts  =  max(max(Tc) ,  tnax(Td)) ;  if  maxpts  <  0  .maxpts  =  0 ;  end 
if  minpts  =  =0  &  maxpts  =  =0,  maxpts  =5,  end 
u  =  [0  pts  1.2*minpts  1.2*maxpta]; 
axis(u); 

subplot<211).plot([0:pts-l],coeffc(index(xl).l:pls).'*') 
titled'Multiphase  Approximation Coeff  at  Resolution  Level  -'.num2str(xl)]) 
axis(u); 

Bubplou;212),plot([0:pts-l],coeffd(index(xl),l:pts),'••) 
title(' Multiphase  Detail  Coefficients') 

xlabel(['n  sampling  points  n*',num2slr<l/f8amp), 'seconds']); 

pause 
strplt 
clg;clc 
axis; 

xl=xl  +  l; 
end 
end 


%  Multi-Phase  Energy  Check 

enrg  1  dmp 

% 

%  Phase  Plane  Determination 

c  =  coeffc.*2; 

d  =  coeffd.*2; 

[i  j]=size(coeffc); 

x  =  0.LL-l; 

yc  =  0:i-l; 

yd=l;i-l; 

cc  =  c(  1 : -lowest+ 1 , 1 :  LL) ; 

clc 

mesh(cc); 

title('Multiphase  Energy  Display  of  the  Approximation') 

pause 

strplt 
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clg 

N  =  10; 

conlour(cc,10.x,yc); 

tilleC Multiphase  Contour  Display  of  the  Approximation'); 

ylabcl('-n  resolution  level'); 

xlabel('10  contour  levels') 

pause 

slrplt 

cl8 

ql  =  inpul('Doyou  desire  to  zoom  in/out  on  the  display  (Y/N)?  W); 

while  ql  =  ='Y'  |  ql   ==  'y', 

q2  =  input('Doyou  want  to  see  the  plots  again  (Y/N)?',V); 

if  q2=  =  'Y'  !  q2  ==  'y', 

mesh(cc) 

tiUe(' Previous  Multiphase  Energy  Display  of  the  Approximation'); 

pause 

contour(cc,N,x,yc) 

ulleC Previous  Multiphase  Contour  Display  of  the  Approximation'); 

ylabel('-n  resolution  level'); 

xlabel(['Numberof  contours:  ',num2str(N)]) 

pause 
end 

clc 

disp(['      ']') 

disp(['x  range:  0  to  ',num2str(j-l)]) 

disp(['y  range:  0  to  -',num2su"(i-l)]);disp(('      ']') 

xl  =input(' Enter  the  minimum  x-value:')  +  l; 

x2  =  input(' Enter  the  maxumum  x-value:')  +  l; 

yl  =  input(' Enter  the  minimum  "magnitude"  y-value  (least  negative):')+  1; 

y2  =  input('Enterthe  maxumum  'magnitude"  y-value  (most  negalive):')  +  l; 

yl  =abs(yl);y2  =  abs(y2); 

y  =  l:l:i;y  =  fliplr(y); 

yll=y(y2); 

y22=y(yl); 

clg 

cc=c(yll:y22,xl:x2); 

mesh(cc) 

litle('Zoomed  Multiphase  Energy  Display  of  the  Approximation'); 

pause 

slrplt 

clg 

N  =  10; 

Nl=10; 

while  NK999 

N  =  N1; 

conlour(cc,Nl,xl:x2,yl-l:y2-l) 

title('Zoomed  Multiphase  Contour  Display  of  the  Approximation'); 

ylabel('-n  resolution  level') 

xlabel(['Numberof  contours:  ',num2str(Nl)]) 

pause 

slrplt 

N 1  =  input( '  Enter  the  number  of  contour  levels  (  <  RET  >  to  continue)  ') ; 
end 
ql  =  input('Zoomin  or  out  Further  (Y/N)?  W); 
end 
clc 

%    Now  the  Detail  Signal 

dd=d(l:-lowe«t,l:LL); 

mesh(dd) 

title(' Multiphase  Energy  Display  of  the  Detail') 

pause 

slrplt 

clg 

contour(dd,10,x,yd) 

utle('Muluphase  Contour  Display  of  the  Detail'); 
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ylabel('-n  resolution  level'); 

xlabel('10  contour  levels') 

pause 

strplt 

clg 

clc 

ql  =  input( ' Do y ou  desire  lo  zoom  in/oul  on  the  display  (Y/N)?','s'); 

while  ql  =  =  'Y'  |  ql  =  =  'y', 

q2  =  input('Doyou  want  to  see  the  plots  again  (Y/N)?',V); 
ifq2=  =  'Y'  |  q2  =='y'. 

mesh(dd);title(' Previous  Multiphase  Energy  Display  of  the  Detail'); 

pause 

conlour(dd,N,x,yd),utle('PreviousMuluphase  Contour  Display  of  the  Detail'); 

ylabel('-n  resolution  level'); 

xlabel([' Number  of  contours:  ',num2str(N)]) 

pause 

end 

disp(['x  range:  0  to  ',num2str(j-l)]) 
disp(['y  range:  -1  to  -',num2str(i-l)]);disp([       ]') 
xl  =input(' Enter  the  minimum  x-value:')  +  1; 
x2  =  input(' Enter  the  maxumum  x-value:')+  1; 

yl  =input(' Enter  the  minimum  "magnitude"  y-value  (least  negative):'); 
if  yl  =  =  0,  yl  =  1;  end; 

y2  =  input(' Enter  the  maxumum  "magnitude'  y-value  (most  negative):'); 
y  1  =  abs(y  1 )  ;y  2  =  abs(y2) ; 
y  =  l:l:i-l;y  =  fliplr(y); 
yll=y(y2); 
y22=y(yl); 
clg 

dd=d(yll:y22,xl:x2); 
mesh(dd) 

title{' Zoomed  Multiphase  Energy  Display  of  the  Detail'); 
pause 
strplt 
clg 
while  N  <  999 

contour(dd,N,xr.x2,yl:y2) 

title('Zoomed  Multiphase  Contour  Display  of  the  Detail'); 

ylabel('-n  resolution  level') 

xlabel(['Number of  contours:  ',num2slr(N)]) 

pause 

•trplt 

N  =  input(' Enter  the  number  of  contour  levels  (<  RET  >  to  end):  '); 
end 
ql  =input('Zoomin  or  out  further  (Y/N)?  W); 
end 


%  Multiple  Phase  Reconstruction 


mprecon 
% 


%       The  End 
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%. .......................... 

%  enrgldmp.m         Determining  the  Energy  in  each  resolution  level 

%  Multi-Phase  One  Dimensional  Case 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%            Phone:     (408)  646-3044/2772 
%.................................... .., 

%  For  the  Multi-phase  case,  the  energy  is  normalized  to  resolution 

%  level  0.    All  the  energies  of  the  "c's"  and  the  "dV  of  a 

%  particular  level  should  add  up  to  the  energy  of  the  *c's"  of  the 

%  next  higher  level.  The  factor  of  two  is  accounted  for  in  the  Mulli- 

%  Phase  case 

% 

%        enrgc  --  energy  array  of  the  "c's" 

%        enrgd  -  energy  array  of  the  "d's" 

%        norm   -  energy  of  the  data  "c0" 

%        esum   -  energy  sum  array  (enrgc + enrgd) 

% 

%  The  calling  routine  is  "mp.m"  . 

%  wavld.m  and  slrplLm  are  also  necessary  m-files 

% 

enrgc  =  zeros(l, -lowest); 
enrgd  =  zeros(  1 ,  -lowest) ; 
a  =  0; 

norm  =  sum(c0  .*2); 
for  row  =  -lowest:-l:l 

a  =  a+l; 

enrgc(row)  =  sum(coeffc(row,:).*2)*2*(-a); 

enrgd(row)  =  sum(coeffd(row,:).*2)*2*(-a); 
end 

enrgc  =  enrgc/norm; 
enrgc  =  [enrgc  1]; 
enrgd  =  enrgd/norm; 
clg 

axis([lowest-2  1  0  1.2]); 
subplot(211),bar([lowest:0], [enrgc]) 

title(['Normalized  energy  of  the  *c'  coefficients    Multi-Phase  Case']) 
xlabeK'Resolution  level') 
axis([lowest-2  1  0  1.2]); 
subplot(2 1 2),bar([lowe«t:  - 1  ] ,  [enrgd]) 
title('Normalized  energy  of  the  "d"  coefficients') 
xlabeK'Resolution  level') 
pause 
strpll 


%  Energy  check  of  the  previous  plot 

esum  =  enrgc  +  [enrgd  0] ; 

clc 

dd  =  ['         Verification  of  the  energy  in  each  resolution  level 


'Level        Energy  in 

Energy  in  C 

Difference 

level  m-1 

level  m 

• 

c+d 

--'1- 

disp(dd) 

J- 

j  = -lowest; 

for  i  =  l:-lowest 

a  =  csum(i)  ;b  =  enrgc(  1 , 

•  +1); 

c  =  l-j;d=abs(a 

-b); 

fprintf('%2.0f       *17.15f 

%17.15f      ' 

.c.a 

.b) 
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fprintfC%17.15f\n\d) 

j=j-i; 

end 

pause 

clc 

clg 

axis; 


%  mpplane.m     develop  Ihe  phase  plane  diagram 

% 

%  By:  LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

%  Phase  Plane  Diagram  for  the  One  Dimensional  Single  Phase  case 

%  This  routine  is  called  by  the  parent  routine  "wavld.m".    All  the 

%  variables  defined  in  "wavld.m"  are  used  here.  New  variables: 

% 

%      c  ( 1 -lowest)  X  LL  "c"  coefficient  matrix 

%       d  (-lowest)  X  LL  "d"  coefficient  matrix 

%       N,N1  —  variables  for  the  number  of  contour  levels 

% 

%  The  main  subroutine  used  is  ptzero.m  funtion  m-file 

% 

clc 
clg 

%    Set  up  the  c  and  d  matrices  from  the  "c*"  and  "d*"  variables 

c=zeros(-lowest+  l.LL); 
d  =  zeros(-  lowest,  LL) ; 
c(l,l:length(cO))=cO; 
clear  cwork  dwork 
for  lvl  =  l:-lowest 

cwork  =  evaI(['ptzero(rmpc',num28tr(lvl),',2*(lvl)-l)']); 

dwork  =  eval(['ptzero(rmpd',num28U-(lvl),',2*(lvl)-l)']); 

c(lvl  + 1 , 1 :  length(cwork))  =  cwork; 

d(lvl,l:length(dwork))=dwork; 
end 

c  =  flipud(c.*2); 
d  =  flipud(d.*2); 


[i  j]=size<c); 

x  =  0:j-l;yc=0:i-l;yd  =  l:i-l; 

mesh(c); 

title<[',MP*  Energy  Display  of  the  Approximation    Phase:  '.. 

setstx<fliplr(phsvct+48))]); 
pause 
N  =  10; 

conlour(c,10,x,yc); 

title(['*MP*  Contour  Display  of  the  Approximation    Phase: 
setstr(fliplr(phsvct-l-48))]); 

ylabel('-n  resolution  level'); 

xlabel('10  contour  levels') 
pause 
clg 
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ql  =  input( ' Do y ou  desire  to  zoom  in  on  the  display  (Y/N)?',V); 
while  ql  =  ='Y"  |  ql  =  =  'y', 
q2  =  input('Doyou  want  to  see  the  plots  again  (Y/N)?'.V); 
ifq2=='Y'  I  q2  ==  y, 

mesh(c); 

title('*MP*  Previous  Energy  Display  of  the  Approximation'); 

pause 

conlour(c,N,x,yc); 

title('*MP*  Previous  Contour  Display  of  the  Approximation'); 

ylabel('-n  resolution  level'); 

xlabel(['Numberof  countours:  ',num2str(N)]) 

pause 
end 
clc 

disp(['     ']') 

disp(['x  range:  0  to  \num28tnj-l)]) 
disp(['y  range:  0  to  -',num2str(i-l)]);disp([       ]') 
xl  =  input( 'Enter  the  minimum  x-value:')  +  l; 
x2  =  input('Enterlhe  maxumum  x-value:')  +  1; 

yl  =  input(' Enter  the  minimum  'magnitude'  y-value  (least  negative): ')+  1; 
y2  =  input(' Enter  the  maxumum  'magnitude*  y-value  (most  negative): ')  +  l; 
yl=abs(yl);y2  =  abs(y2) ; 
y  =  l:l:i;y  =  fliplr(y); 

yll=y(y2); 

y22=y(yl); 

clg 

mesh(c(yll:y22,xl:x2)) 

titie(['*MP*  Zoomed  Energy  Display  of  the  Approximation     Phase:  ',... 

setstr(fiiplr(phsvct + 48))]) ; 
pause 
strplt 
clg 

Nl  =  10; 
while  NK999 
N  =  N1; 

contouKc(yll:y22,xl:x2),Nl,xl:x2,yl-l:y2-l) 
title(['*MP*  Zoomed  Contour  Display  of  the  Approximation    Phase:  ',.. 

setstr(fliplr<ph8vct+48))]); 
ylabel('-n  resolution  level') 
xlabel(['Numberof  contours:  ',num2str(Nl)]) 
pause 
strplt 

Nl  =input(' Enter  the  number  of  contour  levels  (999  to  continue)  '); 
end 

ql  =input('Zoomin  or  out  Further  (Y/N)?  \V); 
end 
clc 

%    Now  the  Detail  Signal 

mesh(d); 

title([',MP*  Energy  Display  of  the  Detail      Phase:  ',... 

setstr(fliplr(ph8vct+48))]); 
pause 
clg 

contoured ,  1 0,  x  ,yd) ; 

title(['*MP*  Contour  Display  of  the  Detail     Phase:',... 

setstr(fliplr(phsvct+48))]); 
ylabel('-n  resolution  level'); 
xlabel('10  contour  levels') 
pause 
clg 

clc 

ql  =input('Doyou  desire  to  zoom  in  on  the  display  (Y/N)?',V); 
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while  ql  =  =  'Y'  |  ql  =  =  'y', 

q2  =  input('Doyou  want  to  see  the  plots  again  (Y/N)?'.V); 
ifq2=  =  'Y'  |  q2  =  =  y, 

mesh(d);title('Previous  Energy  Display  of  the  Detail*); 

pause 

contoured. N,x,yd);title('PreviousContour  Display  of  the  Detail'); 

ylabel('-n  resolution  level'); 

xlabel(['Numberof  contours:  ',num2str(N)]) 

pause 
end 

disp(['x  range:  0  to  ',num2str(j-l)]) 
disp(['y  range:  -1  to  -',num2slr(i-l)]);disp([       ]') 
xl  =inpul('  Enter  the  minimum  x-value:')  +  l; 
x2  =  input( 'Enter  the  maxumum  x-value:')  +  1; 

yl  =  input( 'Enter  the  minimum  'magnitude*  y-value  (least  negative):'); 
if  yl  ==  0,  yl  =  l;  end; 

y  2  =  inpul( '  Enter  the  maxumum  "magnitude"  y-value  (most  negative):'); 
yl  =abs(yl);y2  =  abs(y2); 
y  =  l:l:i-l;y  =  fliplr(y); 
yll=y(y2); 
y22=y(yl); 
clg 

mesh(d(yll:y22,xl:x2)) 
title([  ""MP*  Zoomed  Energy  Display  of  the  Detail     Phase:  ',... 

setstr(fliplr(phsvct +48))]); 
pause 
strplt 
clg 
while  N  <  999 

contour(d(yll:y22,xl:x2),N,xl:x2,yl:y2) 

tiUe([''*'vIP*  Zoomed  Contour  Display  of  the  Detail   Phase:  ',... 
setstr(fliplr(phsvct+48))]); 

ylabel('-n  resolution  level') 

xlabel(['Number  of  contours:  ',num2str(N)]) 

pause 

strplt 

N  =  input( 'Enter  the  number  of  contour  levels  (999  to  end):  '); 
end 

ql  =input('Zoomin  or  out  further  (Y/N)?  W); 
end 


function  x  =  z0pad( Data, h, res lvl) 

*•••• * • 

%  zOpad.m      •*•  For  the  Multi-phase  decomposition  routine  *** 

%  One  Dimensional  Case 

% 

%  By:       LT   J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

% 

%  This  routine  makes  a  row  vector  of  the  lower  coefficients  of 

%  a  particular  resolution  level.    Use  for  only  the  Multi-phase 

%  decomposition  One- Dimensional  case  ! ! 

% 

%  Data  is  the  data  vector 

%  h  is  either  the  h  coeff  or  the  g  coeff  vector 

%  reslvl  is  the  resolution  level 
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m  =  reslvl; 
L  =  length(h); 
J  =  L; 
x  =  0; 
fori  =  l:L 

x  =  x  +  hQ,[zeros(l.(i-l)»2i(ro-l))Datazeros(l.(j-l)*2*(m-l))]; 

end 


function  binary  =  num2bin(number,lengthl) 

% 

%  NUM2BIN.M  Number  to  Binary  conversion 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

% 

% 

%  This  routine  conducts  a  conversion  of  a  base  10  integer  to  a  base  2 

%  data  array  i.e.    [10  10  11] 

% 

%  The  following  variables  are  further  defined 

% 

%        binary      -  binary  output  data  array 

%        lengthl     -  input  desired  length  of  the  display 

%        number      —  base  10  integer 

% 

%  The  following  m-files  are  necessary  for  proper  operation: 

% 

%        wavld.m     --  main  program  (great  grandparent  routine) 

%        mp.m         —  grandparent  routine 

%        mprecon.m  —  multiphase  reconstuction  calling  routine 


binary  =  fj; 
while  number  >  .5 
number = number/2 ; 
if  number- fix( number)  =  =  .5 
number  =  number- .  5 ; 
binary  =  [1  binary]; 
else 

binary  =  [0  binary] ; 
end 
end 

if  length(binary)  <  lengthl 

binary  =  [zeros(  1  ,ab»(  length  1  -length(binary)))binary] ; 
end 
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function  num  =  bin2num(bin) 

* 

%  BIN2NUM.M  Binary  to  number  conversion 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps  navy.mil 

%  Phone:     (408)  646-3044/2772 

% 

% 

%  This  routine  conducts  a  conversion  of  a  base  2  data  array 

%  i.e.    [1  0  1  0  1  1]  to  a  base  10  integer 

% 

%  The  following  variables  are  further  defined 

% 

%        bin      -  binary  output  data  array 

%        num      --  base  10  integer 

% 

%  The  following  m-files  are  necessary  for  proper  operation: 

% 

%        wavld.m     --  main  program  (great  grandparent  routine) 

%        mp.m         —  grandparent  routine 

%        mprecon.m—  multiphase  reconstuction  calling  routine 

% 


num  =  0; 
a  =  length(bin); 
power=a; 
for  ii  =  1 :  a 

num =num  +  bin(ii)*2  "(power-  ii); 
end 


91 


B.      TWO-DIMENSIONAL  ROUTINES 


%. ...................................... ............ 

%  WAV2D.M  15  SEP  92 

%        Wavelet  decomposition  algorithm  for  the  two-dimensional  case 

%        for  either  the  single  selectable  (of  4  possible)  or  the 

%        multiple  phase  case. 

% 

%  By:       LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

or........ •»....••*.**.•***».•......•****.*.**.*...*..*....*.*.*..*.. 

%  This  routine  does  a  two-dimensional  dicrete  wavelet  decomposition 

%  of  a  two-dimensional  data  array.    The  following  routines  are  also 

%  necessary  for  the  proper  operation  of  the  algorithm: 

% 

%  chkputer.m  —   checks  for  the  minimum  hardware 

%  mp2d.m         —   2-D  multiple  phase  decomposition 

%  ptzro2d.m    —   2-D  row  &  col  zero  padding 

%  slrplt.m      -   metafiles  of  selected  plots 

%  hcoefin       -   routine  for  selected  h  coefficients 

%  phas  ??  .m   -    routines  for  phase  00,01,10,1 1  decomp 

%  enrg2d.m      -    energy  check  for  the  2-D  case 

%  recon2d.m     -    2-D  reconstruction 

% 

%    Variables: 

% 

%      pltcnt       -     plot  counter  for  hard  copy  displays 

%      Imagesu"      -     string  variable  of  the  input  data 

%       cO  -     resolution  level  0  data  input 

%       cl  lo  cN      -     corresponds  to  the  c  coeff  level  -1  to  -N 

%      dlltodlN   -     corresponds  to  the  dl  coeffs 

%      d21  to  d2N   -     corresponds  to  the  d2  coeffs 

%      d31  to  d3N   -     corresponds  to  the  d3  coeffs 

%      Hh,Hv,Gh,Gv  -     h  and  g  coeffs  in  the  horiz  &  vert  directions 

%      HhHv.HhGv, 

%       GhHv.GhGv    -     2-D  decomposition  masks 

%       lowest         -     lowest  resolution  level  (a  neg  number) 

%      phsvct?       -     phase  vector  for  either  the  h  or  v  direction 


%  initialization  of  the  counters  and  plots 


clg 

axisCnormal'); 
hold  off; 
pltcnt  =  0; 
chkputer 

clc 


% 

%    2-D  data  input  routine 

Imagestr  =  input( ' Enter  the  name  of  the  image  data:  '  ,'t 
cO  =  eval(Imagestr); 
[il  jl]=size(c0); 

if         jl<  =1  |  il<=l 

disp('Bad  Input  Data,  Restart  the  Program') 
return 
end 


92 


clc 

% 

% 

%  2-D  wavelet  coefficient  input  routine 

bcoefin 

clg 

% 

% 

%  Branch  for  the  2-D  multiphase  case 

disp(['  •]•) 

q  =  input('Do  you  desire  the  2-D  multiphase  decomposition  only?  ',' 

ifq==y  I  q=='Y' 

mp2d 

return 
end 

% 

% 

%  Initialization  for  the  decomposition  routine 

%       of  the  single  selectable  phase 

% 

%  initializing  Cout  for  the  decomposition  routine.  Coul  along 

%  with  Dlout,  D2out  and  D3out  are  intermediate  working 

%  variables  for  the  determining  the  coefficients  for  a 

%  particular  resolution  level  and  phase.  "phas??.m"  are 

%  one  of  for  possible  phase  decomposition  routines  for  the 

%  two  dimensional  case. 

Cout  =  cO; 

subplot(211).me8h(cO),tiUe('3-DPlot  of  the  Data  Array') 

subplot(212),contour(cO),titie('ContourDi8play  of  the  Data  Array') 

pause 

strplt 

clg 

clc 

disp(f  ']') 

q  =  ['Enter  the  number  of  resolution  levels  desired' 

'for  decomposition  of  the  image:  ']; 

disp(q) 

lowest  =  input(' :  >    '); 
clc 
lowest  =  -abs(  lowest) ; 

phsvcth  =  []; 
phsvctv  =  fj ; 

%  The  actual  decomposition  for  the  single 
%        selectable  phase  case 

for  lvl  =  - 1 :  - 1 :  lowest 

clc 

q  =  ['Enter  the  desired  phase  decomposition:    ' 


'A)   Horizontal  Phase  0 

Vertical  Phase  0 

'B)   Horizontal  Phase  0, 

Vertical  Phase  1 

'C)    Horizontal  Phase  1 , 

Vertical  Phase  0 

'D)    Horizontal  Phase  1 

Vertical  Phase  1 

disp(['          ']') 

disp(q) 

phase  =  input( '  :  \V); 

if       phase  =  = 'A' |  phase  == 'a' 
phsvcth  =  [0  phsvcth] ; 
phsvctv  =  [0  phsvctv] ; 
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phasOO 
elseif  phase  =  =  'B'  |  phase  =  =  'b' 

phsveth  =  [0  phsvclh] ; 

phsvctv  =  [l  phsvctv]; 

phasOl 
elseif  phase  =  =  'C  |  phase  =  =  'c' 

phsvclh  =  [l  phsvclh]; 

phsvctv  =  [0  phsvctv] ; 

phaslO 
elseif  phase  =  =  'D'  |  phase  =  =  'd" 

phsveth  =  [1  phsveth]; 

phsvctv  =  [1  phsvctv]; 

p  has  11 
else 

dispC  Error') 

return 
end 

end 

%  clean  up  for  better  memory  utilization 

%  This  is  primarily  for  386  MATLAB 

clear  wrkspc  Dlout  D2out  D3out  Cout  wrk  a  b  q  qq  qqq  ql 

pack 


%  Energy  check  for  the  2-D  single  selectable  phase  case 

enrg2d 

%  clean  up  for  better  memory  utilization 

%  This  is  primarily  for  386  MATLAB 

clear  E  k  lvl  a  b  c  j 

pack 

% 

% 

%  Display  of  the  phase  sequence  in  reverse  bit  order  as  per  the  notes 

clc 

disp(['        •]') 

dispCOrder  of  the  Selected  Phases  from  the  highest  to  the  lowest  level') 

disp('  •) 

disp('Horizontal:  Reading  left  to  right  corresponds  to  going  from  the') 

dispC     higher  resolution  to  the  lower  resolution  level     ') 

disp(['    ']') 

disp(fliplr(phsvcth)) 

disp(['        ']') 

dispCVertical:  Reading  left  to  right  corresponds  to  going  from  the') 

dispC     higherer  resolution  to  lower  resolution  level     ') 

disp(['    ']') 

d  isp(  fliplr(phsvctv)) 

disp(['         ']') 

disp('<RET>  to  continue') 

pause 

clc 

% 

% 

%  Conduct  the  recomposition  routine 
recon2d 

% 

clc 

duP(c     •  ]■) 

dispC  ROUTINE  COMPLETE') 

%      TheEndofWAV2D.M 
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%  HCOEFIN.m       2-D  wavelet  coefficient  input 

% 

%  By:       LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School.  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.np8.navy.mil 

%  Phone:     (408)  646-3044/2772 

% • 

%  This  routine  obtains  the  desired  wavelet  h  coefficients  for  both 
%  the  horizontal  and  vertical  directions.    The  wavelets  are  assumed 
%  separable. 

% 

%  Variables: 

% 

%      Hh,Hv,Gh,Gv  -     h  and  g  coeff  s  in  the  horiz  &  vert  directions 

%      HhHv.HhGv, 

%       GhHv.GhGv    -     2-D  decomposition  masks 

%      wavelet?      -     string  variable  in  the  h  or  v  directions 

% 

* 

%  Obtaining  the  horizontal  h  coefficients 
disp(['  •]•) 

ql  =  ['Enter  the  number  for  the  desired  choice        ' 
'for  the  wavelet  in  the  HORIZONTAL  direction: ' 

(1)  Haar  h  coefficients 

(2)  Daubechies  h  coefficients 

(3)  Own  set  of  h  coefficients 

']; 

disp(ql) 
qqq  =  input('  '); 
if         qqq==3, 

Hh  =  input('  Enter  your  own  set  of  H  coefficients'); 

a  =  8um(Hh);b  =  Hh,Hh'; 

while  a~=l  &  b~=.5 

disp('Yourh  coefficients  are  NOT  compactly  supported!') 

Hh  =  input('Re-enterthe  h  coefficients  or  ctrl-C  to  stop'); 

end 

waveleth  =  [  'Own'] ; 
elseif   qqq==2, 

qq  =  input('Howmany  taps  (2-10)?  '); 

Hh  =  daubdata(qq) ; 

waveleth  =  [ 'Daubechies  ',num2str(qq),'-tap']; 
else      Hh  =  [.5.5]; 

waveleth  =  [  'Haar'] ; 

qqq=i; 

end 
clc 
% 

%  Now  for  the  vertical  case 
disp(['         ']') 

q  =  input('Doyou  desire  the  same  wavelet  in  the  vertical  direction?'.  V); 
ifq=  =  'Y'  |  q==  V 
Hv=Hh; 

waveletv  =waveleth; 
else 

disp(['  ']') 

ql  =['Enterthe  number  for  the  desired  choice 
'for  the  wavelet  in  the  VERTICAL  direction:    ' 

'(1)  Haar  h  coefficients 

'(2)  Daubechies  h  coefficients 

'(3)  Own  set  of  h  coefficients 

']; 

disp(ql) 
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qqq  =  input('  '); 
if         qqq==3, 

Hv  =  input('Enteryour  own  set  of  H  coefficients'); 

a  =  sum(Hh);b  =  Hh,Hh'; 

while  a~  =  l  &  b  —  =  .5 

disp('Your  h  coefficients  are  NOT  compactly  supported! ') 
Hh  =  input('Re-entcrthe  h  coefficients  or  clrl-C  to  Blop'); 

end 

waveletv  =  ['Own']; 
elseif   qqq==2, 

qq  =  input( 'How  many  taps  (2-10)?  '); 

Hv  =daubdata(qq); 

waveletv  =  [ ' Daubechies  ' ,num2str(qq) , '-tap'] ; 
else      Hv  =  [.5.5]; 

waveletv  =  ['Haar'] ; 

qqq=i; 

end 
clc 
end 

% 

% 

%         Generating  the  g  coefficients 
Gh  =  fliplr(Hh); 
fori  =  2:2:length(Hh) 

Gh(l,i)=-Gh(l,i); 
end 

Gv  =  fliplr(Hv); 
fori  =  2:2:length(Hv) 
Gv(l.i)  =  -Gv(l,i); 
end 

% 

% 

%         Generating  the  2-D  masks  for  the  decomposition 

Hv  =  Hv'; 

Gv  =  Gv"; 

%the  2*  factor  here  reduces  the  number  of  multiplication  for 

%each  lower  level  coefficients 

HhHv  =  2*Hv*Hh; 

HhGv  =  2*Gv»Hh 

GhHv  =  2'Hv»Gh 

GhGv=2*Gv»Gh 

NHh  =  length(Hh) 

NHv  =  length(Hv) 

%      End 


%  PHAS00.M       2-D  phase  00  decomposition  of  the  data  array 

% 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.npB.navy.mil 

%  Phone:     (408)  646-3044/2772 

%  This  routine  does  a  phase  00  decomposition  given  the  input  array  -  Cout. 

%  horiz:  0          vert:  0 

% 

%  Variables: 

%  Cout,  D lout- D3out  -    working  variables  for  the  coefficients 

%  wrkspc                   -    working  array  for  the  c  coeffs 
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% 

%  Required  m-files: 

%  wav2d.m      -    calling  program 

%  decompll.m  -    plotting  routine  for  the  res  level 

% 

% 

clc 

[rows  cols]  =  size(Coul); 

%check  if  the  row  is  odd,  if  it  is,  pad  with  one  zero 
if  rem( rows/2, floor(  rows/2))  =  =  0 
Cout  =  [Cout;zeros(  1  .cols)] ; 
rows  =  rows  + 1 ; 
end 

%do  the  same  with  the  columns 
if  rem(cols/2,floor(cols/2))  =  =  0 
Cout  =  [Cout  zeros(rows.l)]; 
cols  =  cob  +  1; 
end 

%  generate  the  workspace 
wrkspc  =  [zeros(NHv-l, cols +  2*NHh-3);zeros(  rows, NHh-l)Cout  zeros(rows.NHh-2);. 

zeros(NHv-l.cols  +  2*NHh-3)]; 
clear  Cout  Dlout  D2out  D3out 

%  Now  operate  on  the  wrkspc  with  the  decomposition  masks 
a=0; 

for  rowl  =  l:2:rows  +  NHv-2 
a  =  a+l; 
b=0; 

for  coll  =  l:2:coU  +  NHh-2 
b  =  b+l; 

wrk  =  wrkspc(rowl:rowl+NHv-l, coll  :coll  +  NHh-l); 
Couu>.b)=  sum(sum(HhHv  .•  wrk)); 
Dlout(a,b)  =  8um(8um(HhGv.*  wrk)); 
D2ouu>,b)  =  sum(sum(GhHv.»  wrk)); 
DSouu^a.b^sun^sumCGhGv.*  wrk)); 
end 
end 
%  Store  the  working  variables  with  the  appropriate  name 
eval(['c',num2su-(abs(lvl)),'  =  Cout;']) 
eval([•d^,num28t^(ab8(lvl)),•=Dloul;•]) 
eval([  'd2 '  ,num2sU<ab8(lvl)) , '  =  D2out; ']) 
eval(['d3',num2str<ab8(ivl)),'  =  D3out;']) 
clg 

%  Display  the  variables  C,  Dl,  D2  and  D3  for  the  resolution  level 
dccomplt 


% 

%  PHAS10.M       2-D  phase  01  decomposition  of  the  data  array 

% 

% 

%  By:       LT  J.  E.  Legaapi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.np8.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

or************************************************************************* 

%  This  routine  does  a  phase  10  decomposition  given  the  input  array  -  Cout. 

%       horiz:  1         vert:  0 

% 

%  Variables: 

%  Cout,  Dlout-D3out  -    working  variables  for  the  coefficients 

%  wrkspc  -    working  array  for  the  c  coeffs 
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% 

%  Required  m  tiles 

%  wav2d.m      -    calling  program 

%  decomplt.m  -   plotting  routine  for  the  res  level 

% 

[rows  cols]=size{Cout); 

%check  if  the  row  is  even,  if  it  is.  pad  with  one  zero 
if  rem(  rows/2, floor(rows/2))  =  =  0 
Cout  =  [Cout;zeros(  1  ,cols)] ; 
rows  =  rows  +  l; 
end 

%odd  column  check 
if  rem(coU/2.floor(coU/2))  -  =  0 
Cout  =  [Cout  zeros(rows,  1)] ; 
cols  =  cols+l; 
end 

%  generate  the  workspace 
wrkspc  =  [zeros(NHv-l, cols  +  2'NHh-4);zeros( rows. NHh-2Cout  zeros(rows,NHh-2); 

zeros(NHv-2,cols  +  2,NHh-4)] ; 
clear  Cout  Dlout  D2out  D3out 
a  =  0; 

for  rowl  =  l:2:rows+NHv-2 
a  =  a+l; 
b=0; 

for  coll  =  l:2:cols  +  NHh-2 
b  =  b+l; 

wrk  =  wrkspc(rowl:rowl  +NHv-l,coll:coll  +  NHh-l); 
Cout(a,b)  =  sum(sum(HhHv  .*  wrk)); 
Dlout(a,b)=8um(8um(HhGv.*  wrk)); 
D2out(a,b)=8um(sum(GhHv.*  wrk)); 
D3out(a,b)=sum(sum(GhGv.*  wrk)); 
end 
end 

lvl2  =  abs(lvl); 

evaKrc'.nun^strdvOX^Cout;']) 
eval([  'd  1 '  ,num2su-(lvl2) , '  =  D 1  out; ']) 
eval(['d2\num2str(lvl2),'=D2out;']) 
eval([  'd3 '  ,num28tr<lvl2) , '  =  D3out; ']) 
clg 
decomplt 


%  PHAS01.M       2-D  phase  01  decomposition  of  the  data  array 

% 

% 

%  By:      LT   J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

%. 

%  This  routine  does  a  phase  01  decomposition  given  the  input  array  -  Cout. 

%       horiz:  0         vert:  1 

% 

%  Variables: 

%  Cout,  Dlout-D3out  -    working  variables  for  the  coefficients 

%  wrkspc  -    working  array  for  the  c  coeffs 

% 

%  Required  m-files: 

%  wav2d.m      -    calling  program 

%  decomplt.m-    plotting  routine  for  the  res  level 
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% 

clc 

[rows  cols)  =  size(Cout); 

%check  if  the  row  is  odd,  if  it  is,  pad  with  one  zero 
if  rem(  rows/2,  floor( rows/2))  —  =  0 
Cout  —  [Cout;zeros(  1  ,cols)] ; 
rows  =  rows+l; 
end 

%check  if  the  columns  are  even 
if  rem(cols/2,floor<cols/2))  =  =  0 
Cout  =  [Cout  zero8(rows.  1 )] ; 
cols  =  cols  +  1; 
end 

%  generate  the  workspace 
wrkspc  =  [zeros(NHv-2,cols  +  2*NHh-3);zeros(rows.NHh-l)Cout  zeros(rows,NHh-2);. 

zero8(NHv-2,cols  +  2»NHh-3)] ; 
clear  Cout  Dlout  D2out  D3out 
a  =  0; 

for  rowl  =  1:2: rows +  NHv-2 
a  =  a+l; 
b=0; 

for  coll  =l:2.coU  +  NHh-2 
b  =  b+l; 

wrk  =  wrkspc(rowl:rowl+NHv-l,coll:coll  +  NHh-l); 
Cout(a,b)=  sum(sum(HhHv  .*  wrk)); 
Dloul(a,b)  =  8um(sum(HhGv.*  wrk)); 
D2out(a,b)  =  8um(8um(GhHv.*  wrk)); 
D3out(a,b)  =  8Um(8um(GhGv.*  wrk)); 
end 
end 
lvl2  =  abs(lvl); 

eval(['c',num28titlvl2),'=Cout;']) 
eval([*d  1 '  ,num2str(lvl2) , '  =  D 1  out; ']) 
eval([  'd2 '  ,num28tr(lvl2) , '  =  D2out; ']) 
eval([  'd3 '  ,num28tr(lvl2) , '  =  D3out; "]) 
clg 
decomplt 


% • * ♦ "•• 

%  PHAS11.M       2- D  phase  11  decomposition  of  the  data  array 

% 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

or  a******************************************* ft**************************** 

%  This  routine  does  a  phase  1 1  decomposition  given  the  input  array  -  Cout. 

%      horiz:  1         vert:  1 

% 

%  Variables: 

%  Cout,  Dlout-D3out  -    working  variables  for  the  coefficients 

%  wrkspc  -    working  array  for  the  c  cocffs 

% 

%  Required  m-files: 

%  wav2d.m      -    calling  prog  nun 

%  decomplt.  m  -    plotting  routine  for  the  res  level 

% 

% 

clc 
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[rows  cols]  =size(Cout); 

%check  of  the  row  is  odd,  if  it  is,  pad  with  one  zero 
if  rem(rows/2,floor(rows/2))  -  =  0 

Cout  =  [Cout;zeros(  1  ,cols)] ; 

rows  =  rows  + 1; 
end 

%do  the  same  with  the  columns 
if  rem(cols/2,floor( cols/2))  —  =  0 

Cout  =  [Cout  zeros(rows.l)]; 

cols  =  cols+l; 
end 

%  generate  the  workspace 
wrkspc  =  [zeros(NHv-2, cols +  2*NHh-4);zeros( rows, NHh-2Pout  zeros(rows.NHh-2); 

zeros(NH  v-2 ,  cols  +  2  *  NHh-4)] ; 
clear  Cout  Dlout  D2out  D3out 
a  =  0; 
for  rowl  =  l:2:rows  +  NHv-2 

a  =  a+l; 

b  =  0; 

for  coll  =  l:2:cols  +  NHh-2 

b  =  b+l; 

wrk  =  wrkspc(rowl:rowl  +  NHv-l,coll:coll  +NHh-l); 

Cout(a,b)=  sum(sum(HhHv  .*  wrk)); 

Dloutfa.b^sumCsumCHhGv.*  wrk)); 

D2ouUa.b)  =  8um(8um(GhHv.*  wrk)); 

D3out(a,b)=8um(eum(GhGv.*  wrk)); 

end 
end 

Ivl2=ab8(lvl); 

eval(['c',num28tr(lvl2),'  =Coul;']) 
eval(['dl  ',num2str<lvl2),'  =Dlout;']) 
eval(['d2',num2sti<lvl2),'  =D2out;']) 
eval(['d3'.num2str<lvl2).'  =  D3out;']) 
clg 
decomplt 


%  DECOMPLT. M   displays  the  coefficients  for  a  particular  level 

% 

%  By:       LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School.  Monterey  CA 

%  e-mail:    legaspi@ece.np8.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

% ' 

%  Variables:  as  defined  in  WAV2D.M  and  PHAS??.M 

% 

%  Required  m- files. 

%  WAV2D.M       parent  routine 

%  PHAS??.M      one  of  four  possible  calling  routines 

%  STRPLT.M      metafile  storage  of  selected  plots 

% 

c'g 

8ubplot(22 1 )  ,contour(Cout) ;  xlabel(  'Cm') ; 

title([ 'wavelet  horiz:  '.wavelelh]); 

subplot(222),contour(Dlout);xlabel('Dlm') 

title([ 'wavelet  vert:  \waveletv]) 

subplot(223),contour(D2out);tiUe('D2m') 

xlabel([ 'Horiz phase:  \setstr(fliplr(phsvcth  +  48)),'  Vert  phase:  ',... 

setsu<fiiplr(pn8vctv +48))]) 
subplot(224).contour(D3out);tJUeCD3m') 
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xlabel(['Resolution  Level  \num2sU"(lvl)]) 

pause 

strplt 

clg 

subplot(22 1 )  .mesh(Cout)  ;xlabel(  'Cm') 

tille(['wavelet  horiz:  '.waveleth]); 

8ubplot(222),mesh(Dlout);xlabel('Dlm') 

title(['wavelet  vert:  '.waveletv]) 

subplot(223).me8h(D2out);liUe<'D2m') 

xlabel(['Horiz  phase:  \setstr(fliplr(phsvcth  +  48)),'  Vert  phase: 

setstr(fliplr<phsvctv +48))]) 
subplol(224),mesh(D3out);liUe('D3m') 
xlabel(['Re8olulion  Level  ',num28tr(lvl)]) 
pause 
strplt 

%       END 


or******************************************************************* 

%  enrg2d.m         Determining  the  Energy  in  each  resolution  level 

%  Two  Dimensional  Case 

% 

%  By:       LT   J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.np8.navy.mil 

%  Phone:     (408)  646-3044/2772 

%  For  the  single  phase  case,  the  energy  is  normalized  to  resolution 
%  level  0.    All  the  energies  of  the  "c's"  and  the  "d's"  of  a 
%  particular  level  should  add  up  to  the  energy  of  the  "c's"  of  the 
%  next  higher  level.  The  calling  routine  is  "wav2d.m"  and  "strplt.m" 
%  is  a  necessary  m-file 


%     Variables  (other  than  defined  in  mp2d.m): 

% 

%        E      -    energy  array  for  [c.dl,d2,d3]  for  eaach  level 

%        esum  -    total  energy  of  each  level:  [c(m)  +  dl +d2+d3]  =c(m+ 1) 

% 

%  Initialization  of  the  energy  sum  for  resolution  level  0 

E  =  zeros(-lowest+  1,4); 

E(l.l)  =  »um(sum(c0.*2)); 

lvl  =  -l; 

k=l; 

%  Calculate  the  energies  for  each  resolution  level 

while  lvl>  =  lowest 

E<k+l,l)=8um(8um((eval(['c',num2elr(abs(lvl))])).*2)); 

E(k+l,2)=sum(sum((eval(['dl',num2str(abs(lvl))]))A2)); 

E(k+l,3)  =  sum(sum((eval(['d2',num2str(ab8(lvl))]))'2)); 

E(k  + 1 ,4)  =  sum(8um((eval([  'd3 '  ,num2str<abs(  lvl))])>  *2)) ; 

lvl  =  lvl-l; 

k  =  k+l; 
end 

E  =  E/E(1,1);     %  normalize  the  energy  to  level  0 
E  =  flipud(E),    %  top  row  is  the  lowest  resolution 
%  Display  of  the  energies  by  coefficients 
clg 

k  =  [lowest-l  1  0  1.1]; 
axis(k);subplot(221),bar(lowe8t:  1 :0,E(: .  1)) 
xlabel('c') 
title([ 'wavelet  horiz:  '.waveleth]) 
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axis(k);subplou;222).bar(loweal:l:0.E(:,2)) 

xlabel('dl') 

title([ 'wavelet  vert:  '.waveletv]) 
axis(k);subplot(223),bar(lowe8t:l:0,E(:.3)) 

tiUe<'d2') 

xlahcK 'single  phase  case') 
axis(k);subplol(224),bar(lowe8t:l:0,E(:,4)) 

liUe('d3') 
pause 
strplt 
clg 
clc 

E  =  E'; 

esum  =  suin(E); 
disp(['        ']') 
dd  =  ['         Verification  of  the  energy  in  each  resolution  level 


'Level        Energy  in 

Energy  in  C 

Difference 

m           level  m-1 

level  m 

c  +  dl+d2+d3 

>l. 

disp(dd) 

J> 

j  =  -lowest; 

for  i=l:-lowest 

a  =  esum(i);b  =  E(l,i+l);c  = 

-j  +  l;d  =  abs(a-b); 

fprintf('%2.0f       %17.15f 

*17.15f      ',c, 

a.b) 

fprintf('%17.15f\n',d) 

i-j-i: 

end 

pause 

clc 

%                    End 

%  recon2d.m        2-D  reconstruction  algorithm  for  the  single 

%  selectable  phase  case 

% 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  Iam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

% 

%      Variables  (other  than  defined  in  wav2d.m): 

% 

%        HvHh.HvGh  -    reconstrution  masks 

%        GvHh.GvGh 

%        whichl        -    phase  decomposition  for  the  particular  level 

%        cwork         -    reconstructed  c  coefficient  array  for  a  level 

% 

%       Required  routines: 

%  wav2d.m 

%  strplt.  m 

%  mmlt??.m   -??=  00,01,10.11  for  the 

%  mask  multiplication 

% 

% 

clc 

disp(['  ']') 

q  =  input('Doyou  desire  to  do  the  reconstruction  (Y/N)?  W); 
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ifq==V  |q=  =  'Y' 

%  The  following  four  lines  generate  the  reconslruclion  masks 

HvHh  =  fliplr(flipud(HhHv)); 

HvGh  =  fliplr(nipud(GhHv)); 

GvHh  =  fliplr(flipud(HhGv)); 

GvGh  =  niplr(nipud(GhGv)) ; 

disp(["   •]') 

q  =  inpul('Doyou  desire  to  sec  the  masks?  ','s'); 

ifq=='Y'  |  q==V 

%  We  will  border  the  arrays  with  zeros  for  a  better  mesh  plot 
[aabb]=size(HhHv); 
dispC  ') 

disp('HhHv  decomposition") 
disp(HhHv) 
subplot(121) 

axis(' square') 

mesh([zeros(  1  ,bb  +  2) ;zeros(aa,  1  )HhHv  zeros(aa.  1 ) ;zeros(  1  ,bb  +  2)]) 

title('HhHv  decomposition  mask') 
dispC  ') 

disp('HvHh  reconstruction") 
disp(HvHh) 
subplot(122) 

axis( 'square') 

mesh([zeros(l,bb  +  2);zeros(aa,l)HvHh  zeros(aa,l);zeros(l,bb  +  2)]) 

title('HvHh  reconstruction  mask ') 
pause 
strplt 
clg 

dispC  ') 

disp('HhGv  decomposition') 
disp(HhGv) 
subplot(121) 

axis('square') 

mesh([zeros(l,bb  +  2);zeros(aa,l)HhGv  zeros(aa.l);zeros(l.bb  +  2)]) 

U'Ue('HhGv  decomposition  mask') 
dispC  ') 

disp('GvHh   reconstruction') 
disp(GvHh) 
subplot(122) 

axis(  'square') 

mesh([zeros(  1  ,bb  +  2);zeros(aa,  l)GvHh  zeros(aa,  1  );zeros(  1  ,bb  +  2)]) 

tit]e('GvHh  reconstruction  mask') 
pause 
strplt 
clg 

dispC  ') 

dispCGhHv  decomposition") 
disp(GhHv) 
subplot<121) 

axis(  'square') 

mesh([zeros(  1  ,bb  +  2) ;zeros(aa.  1  )GhH v  zeros(aa,  1 ) ;zeros(  1  ,bb  +  2)]) 

ut)e('GhHv  decomposition  mask') 
dispC  ') 

disp('HvGh  reconstruction') 
disp(HvGh) 
subplol<122) 

axis(  'square') 

mesh([zeros(l,bb  +  2);zeros(aa,l)HvGh  zeros(aa,l);zeros(l,bb  +  2)]) 

title('HvGh  reconstruction  mask") 
pause 
strplt 
clg 
dispC  ") 

disp('GhGv  decomposition") 
disp(GhGv) 
subplot(121) 

axis('square') 
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mesh([zero8(  1  ,bb  +  2) ;zeros(aa,  l)GhG v  zeros(aa,  1 )  ;zeros(  1  ,bb  +  2)]) 
title('GhGv  decomposition  mask') 
dUp('  ') 

disp('GvGh  reconstruction') 
disp(GvGh) 
subplot<122) 
axis(  'square') 

mesh([zeros(  1  ,bb  +  2) ;zeros(aa,  1  )G vGh  zeros(aa.  1 ) ;zeros(  1  ,bb  +  2)]) 
tiUe('GvGh  reconstruction  mask') 
pause 
strplt 
clg 
end 

clear  aa  bb 
axis('normal') 
clc 

disp(['  ']') 

dispC      NOW  DOING  THE  RECONSTRUCTION') 


%  start  at  the  lowest  level 
%    and  calculate  the  c's  for  the  next  level 
lowest  = -abs(lowest); 
lvl  =  lowest; 

cwork  =  eval(['c',num2str(ab8(lowest))]); 
aa  =  0; 
while  lvl<0 
aa  =  aa  + 1 ; 

[a  b]  =  size(eval([ 'c' ,num2str<ab8(lvl  + 1 ))])) ; 
dl  =eval(['dl'.num28lr(abs(lvi))]); 
d2  =  eval([  'd2 '  ,num2str(abs(lvl))]); 
d3  =  eval([  'd3 '  .num2str(abs(lvl))]); 
which  1  =  [num2str(phsvcth(aa))  ,num2str(phsvcrv(aa))] ; 
eval(['tc  =  mmlt',whichl,'(cwork,HvHh);']) 
eval([  'td  1  =  rnmlt'  .which  1 ,  '(d  1  ,GvHh); ']) 
eval(['td2  =  mmlt',whichl,'(d2,HvGh);']) 
eval(['td3  =mmlt'  .whichl .  '(d3  ,GvGh); ']) 
clear  cwork  dl  d2  d3 

cwork  =  tc(l:a,l:b)+tdl(l:a,l:b)  +  td2(l:a,l:b)+td3(l:a,l:b); 
clear  tc  td  1  td2  td3  a  b  whichl 
clg 

axis(' square') 

%  plot  the  result  of  the  c's  for  comparison 
subplot  121),  contour(  cwork) ; 

title( ' Reconstructed  C  array') 
subplot(  1 22)  ,contour(eval([  'c'  ,num2str(ab8(lvl)- 1 )])) ; 

title([ 'Actual  C  array']) 

xlabel(['Resolution  Level  ',num2str(lvl+  1)]) 
pause 
strplt 
clg 
subplot(  1 2 1 )  ,mesh(cwork) ; 

title( ' Reconstructed  C  array') 
subplot(122),mesh(eval(['c',num28tr(ab8(lvl)-l)])); 

title([ 'Actual  C  array']) 

xlabel( [ ' Resolution  Level  ' , num2str( Ivl  +  1 )]) 
pause 
strplt 
clg 

errr=cwork-eval(['c',nurn28tr(abs(lvl+l))]); 
if  max(max(abs(errr)))  ~  =0 

mesh(abs(enT)) 

titMCAbsolute  error  from  ',num2str(lvl),' to  ',num2str(lvl+  1)]) 

xlabel([ 'Maximum  Error  Value:  ',num2str(max(max(ab8(errr))))]) 

pause 

strplt 
else 

clc 
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clg 

disp(f         ']•) 

dUpC  •*•   NO  ERROR  IN  THE  RECONSTRUCTION   •••') 

disp(['     'J') 

dispf  NOW  CALCULATING  THE  NEXT  RECONSTRUCTION') 

disp(['  FROM  LEVEL  •num28U<lvl  +  l),,TO  ',... 

num2slr(lvl  +  2)]) 
end 

lvl  =  lvl+l; 
end 

axis('normal') 
clc 
end 


function  x  =  mmltOO(array,mask) 

%  MMLTOO.M      mask  multiplication  and  summation  for  the  2-D  phase  00 

%  reconstruction 

% 

%  By:  LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

%    !!!!!!!!  This  routine  is  only  for  PHASE-00  !!!!!!!!! 

%    Required  M  =  files: 

%  wav2d.m     —    grandparent  routine 

%  recon2d.m-   parent  routine 

%  ptzro2d.m  --   puts  n  zeros  rows  &  cols  btwn  coeffs 

% 

array  =  ptzro2d(array ,  1 ) ; 
%  generate  the  workspace 
[NHv  NHh]=size(mask); 
[row  col]  =8ize( array); 

array  =  [array  zeros(row ,NHh-2)] ;cols  =  col  +  NHh-2; 
array  =  [array;  zeros(NHv-2,cols)]; 
a  =  0; 

for  row  1  =  1:  row- 1 
a  =  a+l; 
b  =  0; 

for  coll  =  l:col-l 
b  =  b+l; 

wrk  =  array(rowl:rowl  +  NHv- 1, coll:  coll  +NHh-l); 
x(a,b)  =  8um(sum(mask  *wrk)); 
end 
end 
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function  x  =  mmltlO(array,mask) 

%*•*••»»•.».•••*•»••••*.»•♦•*••♦•♦♦**»•*»•♦*•*•»**♦»*♦♦*..***♦•♦.♦♦*♦*. 

%  MMLT10.M       mask  multiplication  and  summation  for  the  2-D  phase  10 

%  reconstruction 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

^..•..••••••••••. »•••*•*•.»•••»••••••.•».•»*»♦**.•*••♦*♦*••»*••*. •*•••. 

%    !!!!!!!!  This  routine  is  only  for  PHASE-10  !!!!!!!!! 

%    Required  M  =  files: 

%  wav2d.m    -    grandparent  routine 

%  recon2d.m--    parent  routine 

%  ptzro2d.m  --   puts  n  zeros  rows  &  cols  btwn  coeffs 

% 

array  =ptzro2d(array,l); 

%  generate  the  workspace 
[NHvNHh]=size(ma8k); 
[row  col]  =size(  array); 

array  =  [zeros(row,l)  array  zeros(row,NHh-2);zeros(NHv-2,col  +  NHh-l)]; 
a  =  0; 

for  rowl  =  l:row-l 
a  =  a  +  l; 
b  =  0; 

for  coll  =  l:col 
b  =  b+l; 

wrk  =  array(rowl:rowl  +NHv-l,coll:coll  +NHh-l); 
x(a,b)  =  sum(sum(mask.  *wrk)) ; 
end 
end 


function  x  =  mmlt01(array,mask) 

%........ , »..,...., 

%  MMLT01.M      mask  multiplication  and  summation  for  the  2-D  phase  01 

%  reconstruction 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

^**. ...... *..*••. *»••«*«**«*.«.**.«.»*.***...*.*..*««»*«.«*.«.* ........ 

%    !!!!!!!!  This  routine  is  only  for  PHASE-01  !!!!!!!!! 

%    Required  M  =  files: 

%  wav2d.m     --    grandparent  routine 

%  recon2d.m  -   parent  routine 

%  ptzro2d.m—   puts  n  zeros  rows  &  cols  btwn  coeffs 

% 

array  =ptzro2d(array,l); 
%  generate  the  workspace 
[NHvNHh]=size(mask); 
[row  col]  =  size(arTay ) ; 

array  =  [zero8(l,col  +  NHh-2);arrayzero8(row,NHh-2);zero8(NHv-2,col  +  NHh-2)]; 
a  =  0; 
for  rowl  =  l:row 

a  =  a+l; 

b  =  0; 

for  coll  =l:col-l 
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b=b+l; 

wrk  =  array(rowl:rowl  +NHv-l, coll:  coll  +  NHh-l); 
x(a,b)  =  8um(sum(mask .  *  wrk)); 
end 
end 


function  x  =  mmlt  11  (array  .mask) 

%  MMLT11.M       mask  multiplication  and  summation  for  the  2-D  phase  1 1 

%  reconstruction 

% 

%  By:  LT   J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

%    !!!lll!l  This  routine  is  only  for  PHASE- 1 1  !!!!!!!!! 

%    Required  M  =  files: 

%  wav2d.m     —    grandparent  routine 

%  recon2d.m--    parent  routine 

%  ptzro2d.m  -   puts  n  zeros  rows  &  cols  btwn  coeffs 

% 

array  =  ptzro2d(array ,  1 ) ; 

%  generate  the  workspace 
[NHv  NHh]  =size(ma8k); 
[row  col]  =  size(array); 
array  =  [zeros(  1 ,  col)  ;array  ] ;  rows  =  row  + 1 ; 
array  =  [zeros(rows ,  1 )  array] ;  cols  =  col  + 1 ; 
array  =  [array  zeros(rows, NHh-2)] ;cols  =  cols  +  NHh-2 ; 
array  =  [array ;  zeros(NHv-2,cols)]; 
a  =  0; 

for  rowl  =  l:row 
a  =  a+l; 
b  =  0; 

for  coll  =  l:col 
b  =  b+l; 

wrk  =  array(rowl:  rowl  +NHv-l,coll:coll  +NHh-l); 
x(a.b)  =8um(sum(mask.*wrk)); 
end 
end 
function  x  =  ptzro2d(input,n) 


of**************************************************************** 

%PTZR02D.M     puts  "n"  rows  and  columns  of  zeros  after  each 

%  coefficient  in  the  2-D  input  array 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 
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%  Required  routines: 

%  wav2d.m        --  great- grandparent  routine 

%  recon2d.m     -  grand  parent  routine 

%  mmlt??.m      --  parent  routine,  ??  =  00,01, 10,11 

% 

[r  c]  =  size(input); 
xO  =  Q; 
for  j  =  l:c 

x0  =  [x0  input(:,j)  zeros(r.n)]; 
end 

[r  c]=size(x0); 
x  =  D; 
for  j  =  1 :  r 

X  =  [x;x0(j,:);zero8(n,c)]; 
end 


or****************************************************************** 

%  MP2D.M   2-D  Multi-phase  decomposition 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

at****************************************************************** 

%  MP2D  conducts  a  multiple  phase  decomposition  routine  for  the  2-D 

%  case.  The  following  m-files  are  required: 

% 

%        wav2d.m       -    calling  routine 

%        mpdcmp2d.m  -    multiple  phase  decomp  for  1  resolution  level 

%        strpll.m      -   metafile  storage  of  selected  plots 

%        enrg2dmp.m  -    multiple  phase  energy  check 

% 

%    Variables  are  the  same  as  the  calling  program 
% 

% 

% 

%  Initialization  of  the  routine 

clc 

clg 

subplot(2 1 1) .mesh(cO) ,title( "3-DPlot  of  the  Image  Array -) 

subplol(212),contour(cO).tiUe<'ConlourDisplay  of  the  Image') 

pause 

strplt 

clg 

clg 

axis( 'normal') 

hold  off 

dispCf        ']') 

lowest  =  input(' Enter  the  lowest  resolution  level:    '); 

lowest = -abs(  lowest) ; 

m  =  0; 


%    Calculating  the  Coefficients  for  each  resolution  level 

while  m  >  =  lowest  +  1 

eval(['c',num28tr(-m+l),'=mpdcmp2d(c'.num28tr(ab8(m)),',HhHv,m);,]) 
eval([  'd  1 '  ,num28tr(-m  +  1 ) , '  =  mpdcmp2d(c'  ,num2str(abs(m)) , '  ,HhGv  ,m) ; ']) 
eval(['d2',num2str(-m+  1),'  =mpdcmp2d(c',num2str(abs(m)),',GhHv,m);']) 
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eval(['d3',nura2§lr(-m+  l),'=mpdcrnp2d(c',num2str(ab9(m)),',GhGv,m);']) 
%    Display  of  the  coefficients 

subplot(221),eval(['mesh(c',num28tr(abs(m)  +  l ),')']) 
xlabeKCc']) 

title(['wavelet  horiz:  '.waveleth]) 
subplol(222),evald'mesh(dr.num2str(abs(m)+  1),')']) 
xlabel(['dr]) 

titled' wavelet  vert:  ',  waveletv]) 
subplot<223)  .evald  'mesh(d2'  ,num28tr(abs(m)  +  1 ). ') ']) 
title(['d2']) 

\label(  multiphase  case') 
subplot(224)  ,eval([  'mesh(d3 '  ,num2str(abs(m)  +  1 ). ') ']) 
title<[ 'd3']) 

xlabel([' Resolution  level  ',num2su-(m-l)]) 
pause 
strplt 
clg 
subplot(221),eval(['contour(c',num28tr(abs(m)  +  l),')']) 
xlabel(['c']) 

titled 'wavelet  horiz:  '.waveleth]) 
subplot(222)  ,eval([ 'contour(d  1 '  ,nura2str(ab8(m)  + 1 ), ') ']) 
xlabel(['dl']) 

titled  'wavelet  vert:    '  .waveletv]) 
subplot(223).eval(['contour(d2',num2str(ab8(m)  +  l),')']) 
UUe(['d2']) 

xlabel(  'multiphase  case') 
subplou;224).eval(['contour(d3',num28u-(ab8(m)  +  l).,)']) 
UUed'd3']) 

xlabel([ 'Resolution  level  ',num2str(m-l)]) 
pause 
strplt 
clg 
m=m-l; 
end 

% 

% 

%    Energy  check  for  the  2-D  multiphase  case 
enrg2dmp 

% 

*      The  End 


function  x  =  mpdcmp2d(aiTay,mask,m) 
or******************************************************************* 

%  MPDCMP2D.M   2-D  Multi-phase  decomposition  decomposition  function 

%  m-file 

% 

%  By:       LT   J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School.  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.np8.navy.mil 

%  Phone:     (408)  646-3044/2772 

or******************************************************************* 

%  MPDCMP2D  is  the  actual  function  m-file  that  zero  pads  the  input 

%  image  array  given  the  size  of  the  mask  operator  and  the 

%  resolution  level  m.  "m"  is  the  current  resolution  level 

%  and  is  a  negative  number. 

% 

%  mp2d.m        -    calling  routine 

% 

% 
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[NHvNHh]  =  8ize<mask); 

[a  b]  =size(arTay); 

x  =  zeros(a  +  2*(-m)*(NHv-l),b  +  2"(-m)»(NHh-l)); 

fork  =  l:NHh 

for  l  =  l:NHv 

%add  the  rows  of  zeros 

work  =  [zeros(2*(-m)*(l-l),b);array;zeros((NHv.l),2'(-m).b)]; 

[a  b]  =size(work); 

%add  the  coU  of  zeros 

workl  =  [zeros(a,2*(-m)*(k-l))work  zeros(a,(NHh-k),2*(-m))]; 

work2  =  mask(NHv-l  +  1  ,NHh-k  +  l)*workl ; 

x  =  x  +  work2; 
%    clear  work  workl  work2 

end 
end 


%  enrg2dmp.m        Determining  the  energy  in  each  resolution  level 

%  multiphase  case 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

or****************************************************************** 

%  In  this  MULTIPHASE  case,  the  lower  level  will  have 
%  four  possible  phases,  thus  its  energy  will  be  4  times  as  high  as 
%  the  next  higher  level  so  we  will  divide  the  energy  of  the  lower 
%  level  and  the  value  should  equal  the  energy  of  the  c's  of  the 
%  next  higher  level. 

% 

%     Variables  (other  than  defined  in  mp2d.m): 

% 

%        E      -    energy  array  for  [c,dl,d2,d3]  for  eaach  level 

%        esum  -    total  energy  of  each  level:  [c(m)  +  dl +d2+d3]/4  =  c(m+ 1) 

% 

E  =  zeros(-  lowest  +1,4); 

%  Initialization  of  the  energy  for  resolution  level  0 

E(  1 , 1)  =  sum(sum(c0  .  "2)) ; 

lvl  =  -l; 

k  =  l; 

%  Calculate  the  energies  for  each  resolution  level 

while  lvl>  =  lowest 

E(k+l,l)  =  8um(8um((eval(['c,,num2str(ab8(lvl))])).*2*4i(lvl))); 

E(k  +  l,2)  =  sum(sum((eval(['dl,,num2str(abs(lvl))]))*2*4*(lvl))); 

E(k+l,3)  =  sum(sum((eval(['d2,,num2Btr<ab8(lv]))]))"2*4*(lvl))); 

E(k  + 1 ,4)  =  sum(sum((eval([  'd3 '  ,num2str<abs(lvl))]))  i2*4*(lvl))) ; 

lvl  =  lvl-l; 

k=k+l; 
end 

E=E/E(1,1);    %  normalize  the  energy  to  level  0 
E  =  flipud(E);  %  top  row  is  lowest  resolution  level 
%  Display  of  the  energies  by  coefficients 
clg 

k  =  [lowest-l  1  0  1.1]; 
axis(k);subplot(221),bar(lowe8t:l:0,E<:,l)) 
xlabel('c') 
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title(['wavelet  horiz:  \waveleth]) 
axis(k);subplot(222),bar(lowe8t:l:0,E<:,2)) 

xJabel('dl') 

tille(['wavelet  vert:  \waveletv]) 
axis(k);subplot(223),bar(lowest:  1 :0.E(:  ,3)) 

tiUeCd2') 

xlabel('multiphase  case') 
axis(k);subplot(224),bar(lowest:l:0,E<:,4)) 

UUeCdS') 
pause 
strplt 
clg 
clc 

E  =  E'; 

esum  =  sum(E); 
esum  =  sum(E) ; 
disp([*       "J') 
dd  =  ['         Verification  of  the  energy  in  each  resolution  level 


'Level        Energy  in  Energy  in  C  Difference' 

i  level  m-l=  level  m 

c+dl+d2+d3 

']; 

disp(dd) 

j  =  lowest; 

for  i=  l:-lowest 

a  =  esum(i);b  =  E(l,i+l);d  =  abs(a-b);c=j+  1; 

fprintf('%2.0f       %17.15f         %17.15f      ',c,a.b) 

fprintf('%17.15An',d) 

j=j+i; 

end 

pause 

clc 

%  End 


eg************************************************************************* 

%  IDATA.M  Sample  Image  data  for  testing  the  routines 

% 

%  By:       LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

K=menu('Select  image', 'Boxr,'Box2', 'Bo xS'.'BAV'.^BAV'.'gB/W',... 
'X'.'Hexagon'); 


ifK=  =  l 

%  first  image  is  a  box 

im=zeros(100,100); 

im(20:21, 20:80)  =  ones(2, 61) 

im(79:80,20:80)  =  ones(2,61) 

im(22:78,20:21)  =  ones(57,2) 

im(22:78,79:80)  =  ones(57,2) 

elseif  K==2 
im=zeros(100,100); 
im(20,20:79)  =  ones(l,60); 
im(80,20:79)  =  ones(l,60); 
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im(20:79,20)=ones(60,l); 
im(20:79.80)  =  onc8(60.1); 

elseif  K==3 

im  =  zero8(100,100); 

im(l:2,l:61)=ones(2,61); 

im(60:61.1:61)  =  ones(2,61); 

im(3:59,l:2)  =  ones(57,2); 

im(3:59,60:61)  =  ones(57,2); 

eUeif  K==4 

%  second  image  is  a  1  white  and  1  black  rectangle 

im=zeros(100,100); 

im(:,51:100)=ones(100,50); 

elseif  K=  =5 

%  third  image  are  two  black  boxes  and  two  white  boxes 

im=zeros(100,100); 

im(l:50,l:50)  =  ones(50,50); 

im(51 :  100.51 :  100)  =one*(50,50); 

eUeif  K==6 

%  fourth  image  are  8  black  and  8  white  boxes 

im=zeros(100,100); 

im(l:25,:)  =  [ones(25,25)zero8(25,25)one«(25,25)zeros(25.25)]; 

im(26:50,:)  =  [zero8(25,25)ones(25,25)zeros(25,25)ones(25.25)]; 

im(5 1 : 75, :)  =  [ones(25.25)zeros(25,25)  ones(25 ,25)  zeros(25,25)]; 

im(76:100,:)  =  [zero8(25,25)ones(25,25)zero8(25,25)ones(25,25)]; 

eUeif  K==7 

%  fifth  image  are  2  45  degree  diagonals 

im  =  eye<100,100); 

im  =  sign(im  +  flipud(im)); 

eUeif  K  =  =8 

%  sixth  image  is  a  hexagon 

for  row  =  1:20 

n  =  2*row-l; 

im6(row,n:n  +  2)  =  ones(l,3); 
end 

[a  b]  =size(im6); 
row2  =  ceil(sqrt(a*2  +  b*2)); 
im6a  =  zeros(row2,b); 
im6a(:  ,b- 1  :b)  =  ones(row2,2); 
rim  =  [im6;im6a;fliplr(im6)]', 
lim  =  fliplr(rim); 
im6  =  [lim  rim]; 
clear  row  row2  im6a  lim  rim 
[a  b]  =size(im6); 
im6a  =  zeros(100,100); 
aa  =  noor((100-a)/2); 
bb  =  floor((100-b)/2); 
im6a(aa  + 1 :  aa  +  a,bb  + 1 :  bb  +  b)  =  im6; 
im  =  im6a; 

clear  im6  im6a  a  aa  b  bb  K  n 
end 
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C.      ROUTINES  COMMON  TO  BOTH  DIMENSIONS 


%  chkputer.m        checks  for  the  proper  operating  system 

% 

%  By:  LT   J.  E.  Legaspi 

%  Professor  Lam 

%  US  NavaJ  Postgraduate  School,  Monterey  CA 

%  e-mail:    legaspi@ece. nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 


•*•••**••••**•••••••••*•••***•••••*••••**••***•**••*••**••••••«•••• 


% 


%  This  routine  reads  the  variable  "cmptr"  to  see  if  it  can  run  in 

%  the  current  system.    At  this  time,  the  algorithm  will  run  on  a 

%  PC386,  UNIX,  or  SUN  system.    If  the  operating  system  is  a  VAX- 

%  based,  the  routine  will  not  be  as  reliable  since  the  VAX  MATLAB 

%  does  not  recognize  IEEE  arithmetic .    The  routine  also  deletes  all 

%  other  leg*. met  files  prior  to  running. 

% 

[cmptj  mxsz]  =  computer; 

a  =  length(cmptr) ; 

clc; 

if       cmptr<l)  =  ='P'&cmpu-(2)  =  =  'C,&cmptr<3)==,3', 

!del  leg*. met; 
elseif  cmptr(l)  =  =  'S'&cmptr(2)  =  =  'U'&cmptr(3)  =  ='N\ 

!rm  leg*. met; 
elseif  cmptr(a-3)=  =  'U'&cmpu-(a)  =  =  'X', 

!rm  leg*. met; 
else 

clc 

disp(['  •]') 

ql  =['*****Errorin  input******' 
'Hardware  cannot  support   ' 
'the  algorithm... PLS 
'  Restart  the  Program        ']; 

disp(ql) 

return 
end 


% 

%  su-plt.m 

% 

%  Store  Plot  M-file  for  wavelet  decomposition  routine 

% 

%  By:      LT  J.  E.  Legaspi 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    lega8pi@ece.nps.navy.mil 

%  lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044/2772 

% 

%  This  program  stores  the  plot  as  a  metafile  with  the  prefix  "leg" 
%  followed  by  a  number  starting  from  0  on  upwards  and  ending  with 
%  the  extension  ".met"    You  must  be  familiar  with  GPP  to  get  hard- 
%  copies  of  the  plots,  'pltcnl"  is  the  only  variable  necessary  for 
%  the  proper  indexing  of  the  leg*. met  files.    Ensure  that  the 
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%  calling  program  docs  not  redefine  the  variable! 
% 

disp(['    ']•) 

q5  =  inputCStorethe  Plot  (Y/N)?  \V); 

ifq5  ==  T  I  q5==  Y. 

eval(['meta  leg',num2str(pltcnt)]) 

pllcnt  =  pllcnt+l; 

end 


act**************************************************************** 

%  daubdata.m 

% 

%  By:      CAPT/USMC  Kalrabach,  M. 

%  Professor  Lam 

%  US  Naval  Postgraduate  School,  Monterey  CA 

%  e-mail:    lam@ece.nps.navy.mil 

%  Phone:     (408)  646-3044 

% 

% 

function[hcoeff]  =  daubdata(x) 

%    This  function  returns  the  proper  "h"  coefficients  for  the  specified  order 

%    Daubechies  wavelet.  The  input  value  of  "x"  is  the  specified  order. 

if  x  ==  2 

hcoeff  =  [0.482962913145;0.836516303738;0.224143868042;-0. 129409522551] 
elseif  x  =  =  3 
hcoeff  =  [0.332670552950;0.806891509311;0. 4598775021 18;-0. 13501 1020010; 

-0.085441273882;0.035226291882]; 
elseif  x  =  =  4 
hcoeff  =  [0.230377813309;0. 714846570553;0.630880767930;-0.027983769417; 

-0. 18703481 1719;0. 030841381836;0. 03288301 1667;-0. 010597401785]; 
elseif  x  =  =  5 
hcoeff  =  [0.160102397974;0.603829269797;0.724308528438;0. 138428145901; 

-0.242294887066;-0.032244869585;0.077571493840;-0. 006241490213; 

-0.012580751999;0.003335725285]; 
elseif  x  =  =  6 
hcoeff  =  [0. 111540743350;0.494623890398;0. 751 133908021  ;0.3 15250351709 

-0.226264693965;-0.129766867567;0.097501605587;0. 027522865530; 

-0.031582039318;0.000553842201;0.004777257511;-0.001077301085] 
elseif  x  =  =  7 
hcoeff  =  [0.077852054085;0. 3965393 19482;0.729132090846;0.469782287405 

-0.143906003929;-0.224036184994;0.071309219267;0. 080612609151; 

-0. 038029936935;-0. 016574541631  ;0. 012550998556,0.000429577973; 

-0.001801640704;0.000353713800]; 
elseif  x  =  =  8 
hcoeff  =  [0.054415842243;0.312871590914;0.675630736297;0. 585354683654; 

-0.0158291O5256;-0.284015542962;0.0OO472484574;0. 128747426620; 

-0. 017369301002;-0. 044088253931  ;0.013981027917;0. 008746094047; 

-0.004870352993;-0. OOO391740373;0. 00067544  9406;-0. 0001 17476784] 
elseif  x  =  =  9 
hcoeff  =  [0.038077947364;0.243834674613;0. 604823 123690:0.657288078051; 

0.133197385825;-0.293273783279;-0. 096840783223:0. 148540749338; 

0.030725681479;-0. 067632829061  ;0. 0002509471 15;0.022361662124; 

-0.004723204758;-0. OO4281503682;0. 001 847646883:0.000230385764; 

-0.000251963189:0.000039347320]; 
elseif  x  ==  10 
hcoeff  =  [0. 026670057901;0. 188176800078;0.527201188932;0. 688459039454; 

0.281 172343661  ;-0.249846424327;-0.195946274377;0. 127369340336; 

0.093057364604;-0.071394147166;-0. 029457536822;0  033212674059; 

0.003606553567;-0. 010733175483,0. 001395351747;0. 001992405295; 

-0.000685856695;-0.000116466855;0.000093588670;-0.000013264203] 
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else 

hcoeff  =  0; 

break 
end 

hcoeff  =  hcoeff /sqrt(2);  %  normalize  the  "h"  vector 

return 


#!  /bin/csh  -f 

9 

^•••••••••••••••••••••••••••••••♦•••••♦••♦•••••••••••••••••••••••••******** 

§  metal4  C-Shell  Routine 

# 

#  metal4  is  a  command  executable  Tile  for  the  generation  of 

#  the  leg*. met  metafiles  into  postscript  graphic  files  for  the 

#  SPARC  printer  number  14  on  the  fourth  floor  of  Spanagel431. 

#  This  routine  is  system  dependent  and  is  guaranteed  to  work 

#  on  the  UNIX  based  system  of  the  ECE  Department  at  the  Naval 
i  Postgraduate  School. 

# 

Is  -1  leg*. met  >tmmmpl 

echo  "#!  /bin/csh  -f  "  >  tmmmp2 

awk  -F.  '{printf  "gpp  "  $1  "  -dps  -ol  \n"  ;\ 

printf  'lprl4  "  $1  ".ps  \n";\ 

printf 'rm  "  $1  '.ps  \n"}'  tnunmpl  >  >tmmmp2 
chmod  0700  tmmmp2 
tmmmp2 
rm  tmmmpl  tmmmp2 
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